Gérard Eberl
email: gerard.eberl@pasteur.fr

Thomas Pradeu
email: thomas.pradeu@u-bordeaux.fr

Towards a general theory of immunity?

Theories are indispensable to organize immunological data into coherent, explanatory, and predictive frameworks. Here we propose to combine different models to develop a unifying theory of immunity, which situates immunology in the wider context of physiology. We believe that the immune system will be increasingly understood as a central component of a network of partner physiological systems that connect to maintain homeostasis.

Why a theory of immunity?

Since its inception in the late 19th century, immunology has been embedded into medicine, aiming to explain how we do not succumb to infectious microbes.

Immunology is practical: based on our understanding of the immune system, we have developed vaccines and immunotherapies against infection and cancer, and are countering its destructive effects during pathological inflammation. Nonetheless, immunology was awarded several Nobel Prizes for its theories, in 1908 to Metchnikoff and Ehrlich for their complementary frameworks ("cellular" and "humoral"), in 1960 to Burnet for founding concepts on tolerance and the clonal selection theory, in 1984 to Jerne for theories "concerning the specificity in development and control of the immune system". These awards underscore the fundamental challenge in explaining how the immune system recognizes pathogens.

Immunological theories allowed designing experiments to advance our understanding of these mechanisms, and develop tools to translate this understanding into medicine.

Why a new theory?

The self-non-self theory has dominated immunology ever since it was proposed by Burnet in the 1940s, drawing the line between us and infectious microbes by excluding all microbes. The Danger model, developed by Matzinger in the 1990s, proposes however that the immune system is not activated by non-self, but rather by the damage that microbes, or anything else, inflicts to our cells and tissues [START_REF] Matzinger | Tolerance, danger, and the extended family[END_REF].

Nevertheless, microbial products, such as LPS, can activate the immune system without damaging cells and tissues (even though it may be viewed as a proxy of danger). Therefore, Jaeger, Vivier and one of us proposed in 2013 that the immune system is not fundamentally activated by non-self or danger, but rather by a change in normality, termed "discontinuity": the rate of change determines to the amplitude of effector responses [START_REF] Pradeu | The speed of change: towards a discontinuity theory of immunity?[END_REF][START_REF] Pradeu | The discontinuity theory of immunity[END_REF]. Of note, discontinuities come in many types that define the types of immune response, and include anything from incoming microbes, chemicals, damaged extracellular matrix, different modes of cell death or increased expression of self-proteins. This theory generalizes the function of immunity to the maintenance of homeostasis, rather than only to defense against pathogens, a notion discussed previously by and by one of us [START_REF] Eberl | A new vision of immunity: homeostasis of the superorganism[END_REF] in the context of the evolution of microbiota-host interactions. Immunity may be broadened to include tissue, cell and molecular repair, as proposed in 2016 by Kourilsky [START_REF] Kourilsky | The natural defense system and the normative self model[END_REF]. So what is now missing to understand immunity?

The equilibrium model of immunity

These theories explain how the immune system is activated, but do not account for more general levels of immune cross-regulations. Immune responses are induced by discontinuities but, importantly, they come into several mutually inhibitory types [START_REF] Eberl | Immunity by Equilibrium[END_REF].

Type 1 responses are directed against infected or transformed cells, induced by mechanisms that detect intracellular discontinuities, such as the MHC class I antigen processing pathway and cytoplasmic pattern recognition receptors (PRRs). Type 3 responses, in contrast, are induced by discontinuities that affect the extracellular space, such as bacterial or fungal infections, and typically involve phagocytes. In the case of large extracellular threats, such as helminthes, which cannot be eliminated by phagocytes, type 2 responses are induced that use tissue repair processes to block entry of the parasites into tissues. Mucus production and collagen deposition are characteristic of such responses, which may lead to tissue fibrosis if prolonged. A fourth type of response is proposed that prevents microbes from approaching sensitive tissues where they would trigger damaging inflammation, such as the eye or the intestinal lamina propria [START_REF] Matzinger | Tissue-based class control: the other side of tolerance[END_REF].

As we now realize, all these responses are active in a healthy organism, as viruses, bacteria, fungi and parasites are ever present in the microenvironment at mucosal surfaces and within, and tissues have to be constantly repaired. We also know that the different types of immune responses, including regulatory T cells (Tregs), are mutually inhibitory, a cross-regulation that defines a homeostatic equilibrium. An ongoing anti-helminth type 2 response, for example, blocks type 1 anti-viral responses or type 3 anti-bacterial responses [START_REF] Osborne | Coinfection. Virus-helminth coinfection reveals a microbiota-independent mechanism of immunomodulation[END_REF]. Conversely, the antibioticmediated destruction of bacteria leads to decreased type 3 responses, and thus to increased susceptibility to allergy and improved anti-viral responses [START_REF] Baldridge | Commensal microbes and interferon-lambda determine persistence of enteric murine norovirus infection[END_REF].

Immunity, physiology and ecology

The immune system is embedded in an organism that has limited resources [START_REF] Schulenburg | Introduction. Ecological immunology[END_REF]. We all have experienced the exhausting effects of an infection, the immunity to which requires at times massive mobilization of energy. Conversely, stress can divert energy away from the immune system to accomplish an urgent task, such as the fight or flight response to a predator. Thus, the immune system is dependent on the availability of resources stored in cells and specialized tissues, and on resource allocation centrally managed, in many animals, by the brain. The nervous system, in its own way, perceives discontinuities within the organism and in the environment, and coordinates behavioral and physiological processes to maintain homeostasis [START_REF] Kotas | Homeostasis, inflammation, and disease susceptibility[END_REF], and thus, survival of the organism. As emphasized recently by Veiga-Fernandes and Freitas, the immune and the nervous systems are comparable in that both sense perturbations in the organism and in its environment [START_REF] Veiga-Fernandes | The S(c)ensory Immune System Theory[END_REF], and interact intimately through neuroimmune cells "units" [14].

Towards a general and practical theory of immunity

In conclusion, in order to best describe the functioning of the immune system, a general theory of immunity has to take into account at least 3 levels of activation and regulation that orchestrate homeostasis of the organism (Figure 1): 1) the activation of the immune system by discontinuities induced by microbes, tumors, and injuries, 2) the mutual regulation of immune responses to distinct types of discontinuities, and 3) the global regulation of the immune system within the organism. While most of 20 th century immunology has focused on understanding level 1, immunology now integrates microbiota and tissue repair as essential components of its reactivity (level 2). Level 3 is progressively opening up, as physiology and neurosciences develop powerful new experimental tools, while immunologists become aware of the fundamentally transversal nature of their field [15] (Textbox). If the perspective presented here is correct, we can predict that immunology will be increasingly understood as a component of a large network of partner biological domains (in particular endocrinology, neurosciences, the study of metabolism, the study of repairregeneration), in which immunology is likely to play a central role.

Textbox: Assessing the physiological integration of the immune system

Integration of the immune system with the other systems in the organism (see Figure 1, level 3) is necessary to maintain homeostasis. All systems cross-talk to generate energy, growth, movement and defense, while keeping life parameters within physiological levels. We now come to realize that the immune system is key to the function of the other systems and in turn, requires input from these other systems to function. Therefore, a more holistic approach to immunology appears necessary to develop a full understanding of the immune system. Here are a few examples of physiological cross-regulations with the immune system that are coming to the forefront of research, and carry much potential for discovery and progress, both in the lab and in the clinic.

1. Nervous system: the cross-talk between the nervous and immune systems has been investigated for decades, but novel technologies to manipulate both molecules and cells with exquisite precision in both systems, such as optogenetics, allow for detailed mechanistic investigations of nervous and immune networks (Gradinaru, V. et al.,

Science 2009). Recent advances show how immune cells and neurons interact during

health and disease (Kipnis, J., Science 2016), but the regulation of cognition by immunity, and immunity by cognition, remain intriguing cross-talks to explore.

2. Digestive system: the role of the intestinal microbiota in digestion, of its metabolites in the immune system and hepatic functions, the regulation of the microbiota and hepatic functions by the immune system, and, more generally, the cross-talk between metabolism and immunity, rise profound interest in the research community, as perversion of these cross-talks are associated with metabolic diseases such as diabetes ( ). Less clear in many cases is the logic of such regulation at the organismal level, and the effect of the immune system on these processes that control physiological homeostasis during adaptation to the environment.

5. Aging: with time, integration of the immune system may become deregulated as its cross-talk with other systems may suffer from the senescence of stem cells and the consequent loss of regenerative capacity in its own cells and tissues, as well as in cells and tissues of the other systems (Oh, J., Nat Med 2014). Such loss of integration may favor the development of chronic inflammatory pathologies that are characteristic of old age.

14

  Veiga-Fernandes, H. and Mucida, D. (2016) Neuro-Immune Interactions at Barrier Surfaces. Cell 165 (4), 801-811 15 Qiu, Y. et al. (2014) Eosinophils and type 2 cytokine signaling in macrophages orchestrate development of functional beige fat. Cell 157 (6), 1292-1308

  Ghesquiere, B. et al., Nature 2014). Digestion, metabolism and immune reactions are also linked as sources and sink of energy (O'Neill, L.A. et al.,Nat Rev Immunol 2016). In animals, the management of energy at the organismal level requires the nervous system to integrate other energetic priorities, such as movements generated by the skeletal system, itself under the control of perception by the CNS. Thus, exploring energy management requires many systems to be investigated as integrated networks. 3. Skeletal and circulatory systems: bones, muscles and vessels are prone to injury, and must be repaired. Type 2 immune responses co-opt repair mechanisms to fight large parasites or develop allergic reactions (Gause, W.C. et al., Nat Rev Immunol 2013). They also play a central role in tissue repair, the deregulation of which leads to fibrosis and vascular diseases (Gieseck, R.L., et al., Nat Rev Immunol 2017). It remains however unclear to what extent tissue repair processes are dependent on type 2 responses, and to what extent repair processes have shaped pro-allergic and anti-

	helminth responses.

4. Endocrine system: hormones mediating stress, sex, energy partition or circadian rhythms all have effects on immune cells, involving direct effects of nuclear hormone receptors in the regulation of immune genes (Glass, C.K. and Ogawa, S., Nat Rev Immunol 2006


Julien Baptiste

Emmanuel Lepinette
email: emmanuel.lepinette@ceremade.dauphine.fr

Diffusion equations: convergence of the functional scheme derived from the Binomial tree with local volatility for non smooth payoff functions

Keywords: Binomial tree model, European option pricing, Diffusion partial differential equations, Finite difference scheme, Finite element scheme. 2000 MSC: 60G44, G11-G13

Diffusion equations: convergence of the functional scheme derived from the Binomial tree with local volatility for non smooth payoff functions.

Introduction

There are three major methods to simulate a diffusion equation of type

h t (t, x) + x 2 σ 2 (t, x) 2 h xx (t, x) = 0, t ∈ [0, T ), (1.1)
with the boundary condition h(T, x) = g(x). The first one is to use the Monte Carlo methods as the solution h admits a probabilistic representation, see [START_REF] Friedman | Stochastic differential equations and applications[END_REF].

The second one is to numerically compute the solution directly from the PDE.

In particular, the very well known finite difference (FD) technique is based on approximations of the successive derivatives. It requires a discretization, both in time and in the space variable, of some compact subset [0, T ] × [a, b], a ≤ b. Moreover, as we need to fix extra boundary conditions, when x = a or x = b, that are a priori approximations of h, a second type of approximation errors may appear. For instance, if b is large enough, we generally set the condition h(t, b) = lim x→∞ h(t, x) if the latter may be estimated while the condition h(t, 0) = g(0) is chosen if a = 0 under some mild conditions on σ.

There is a lot of articles in the literature focusing on this technique, see for example [START_REF] Brenner | Besov spaces and applications to difference methods for initial value problems[END_REF] and [START_REF] Richtmyer | Difference methods for initial value problems[END_REF]. The finite element (FE) method is more sophisticated, see [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF], but it is based on a discretization of the space variable as well as other techniques as finite volume [START_REF] Eymard | Finite volume methods[END_REF] or spectral methods [START_REF] Quarteroni | An introduction to spectral methods for partial differential equations[END_REF], [START_REF] Bernardi | Spectral methods[END_REF].

The third one is based on binomial trees. Such a scheme is proposed by Milstein [START_REF] Milstein | The Probability Approach to Numerical Solution of Nonlinear Parabolic Equations[END_REF] where a convergence theorem is proven when the terminal condition, i.e. the payoff function in finance, is smooth. In particular, it is supposed that the successive derivatives of the solution h of the P.D.E. are uniformly bounded. In finance, e.g. for the Call payoff function g(x) = (x -K) + , this is not the case. On the contrary the successive derivatives of the P.D.E. solution explode at the horizon date [START_REF] Lépinette | Approximate hedging in a local volatility model with proportional transaction costs[END_REF]. Motivated by pricing in finance, we consider non smooth payoff functions and we study the convergence of this very simple alternative scheme derived on the binomial tree with local volatility financial model, which is only based on a discretization in time of the interval [0, T ] and allows to simulate t → h(t, x) for a fixed x. The advantage of this method is the ease of implementation of the functional scheme which is basically defined. The numerical scheme we study is defined by two functions depending on the number n of discretization dates. We focus here on a particular choice of such functions which leads to a uniform approximation in time of the diffusion equation under some mild conditions. A convergence theorem which generalizes [START_REF] Milstein | The Probability Approach to Numerical Solution of Nonlinear Parabolic Equations[END_REF] to non smooth terminal conditions is proved. We present some examples where the functional scheme appears to outperform the finite difference method and, moreover, it is well adapted to parallel computing. An open problem is to study the suggested functional scheme more generally.

Motivation and interpretation in finance. Let us consider a binomial model in finance with T steps defined as follows: we consider on a stochastic basis (Ω, (F t ) t=0,••• ,T , P) a price process (S t ) t=0,••• ,T where S 0 is given at time 0. Moreover, we suppose that P(S t+1 = S t k u t (S t )|F t ) = P(S t+1 = S t k d t (S t )|F t ) = 1/2 where 0 ≤ k d t < 1 < k u t are two function-valued processes such that

k d + k u = 2.
By the assumptions, we easily deduce that (S t ) t=0,••• ,T is a P-martingale.

Let us define a portfolio process V as an adapted process satisfying ∆V t :=

V t -V t-1 = θ t-1 ∆S t for all t = 1, • • • , T , where θ t is F t -measurable for all t = 0, • • • , T -1.
Moreover, we say that it replicates the contingent claim g(S T ), for some function g, if V T = g(S T ). If g(S T ) is integrable, we easily deduce that V is a martingale and, by induction, we get that V t = C(t, S t ) for some measurable function C. Using the martingale property V t-1 = E(V t |F t ), we finally deduce that

C(t -1, S t-1 ) = 1 2 C(t, S t-1 k d t-1 ) + 1 2 C(t, S t-1 k u t-1 ), t = 1, • • • , T. (1.2)
The natural idea is then to consider this scheme as the number n+1 of dates (T i/n) i=0,••• ,n tends to ∞ when we discretize the continuous-time interval [0, T ]. This leads to the scheme we propose to study in the next section when k d t (x) = 1 -σ(t, x) T /n and k d t (x) = 1 + σ(t, x) T /n. Actually, the asymptotic behaviour of such binomial model is well known, at least for the price process. Indeed, let us consider the continuous-time price process S n of the binomial model defined by

S n t = S n T (i-1)/n if t ∈ [T (i -1)/n, T i/n) where S n
T (i-1)/n is defined from k d and k u when starting from S 0 at t = 0 as explained above. Then, by [START_REF] Prigent | Weak convergence of financial markets[END_REF]Proposition 3.2.1] under mild conditions on σ, the sequence S n weakly converges to the diffusion process S satisfying the stochastic differential equation

dS t = S t σ(t, S t )dW t , (1.3)
where W is a standard Brownian motion. We also deduce the convergence of the price functions (C n ) n≥1 given by (1.2 ) towards the limit price C of the continuous-time model given by (1.3), i.e. the solution to the diffusion equation (1.1). This is indeed the case when σ is constant for a large class of payoff functions g as proven in [2, Chapter 4] where it is also proven that the best rate we can get in general is 1 √ n , see also [START_REF] Walsh | The rate of convergence of the binomial tree scheme[END_REF] and the papers [START_REF] Cox | Option pricing: a simplified approach[END_REF], [START_REF] Hsia | On binomial option pricing[END_REF], [START_REF] Diener | Asymptotics of the price oscillations of a European call option in a tree model[END_REF], and [START_REF] Viitasaari | Rate of convergence for discrete approximation of option prices[END_REF] among others and, finally, [START_REF] Milstein | The Probability Approach to Numerical Solution of Nonlinear Parabolic Equations[END_REF] for smooth terminal conditions and local volatility function. Actually, our main contribution is to show that the result is still valid for local diffusion coefficients even if the solution h of the diffusion equation admits unbounded derivatives, e.g. when the terminal condition is not smooth. Moreover, the numerical experiments we present confirm the accuracy of the functional scheme in comparison to the (FD) and (FE) methods.

Functional scheme for diffusion equations

Let us consider a bounded diffusion function σ : [0, T ]×R + and the associated backward parabolic equation

h t (t, x) + x 2 σ 2 (t, x) 2 h xx (t, x) = 0, t ∈ [0, T ), h(T, x) = g(x), (2.4)
where the terminal condition is defined by a Lipschitz function g with Lipschitz constant L g and h t and h xx are respectively the first and the second derivatives of h with respect to time t and space variable x ∈ R. This equation is very well known in physics but also in mathematical finance for models without friction when the risky asset S is driven by a standard Brownian motion W so that it satisfies the stochastic differential equation dS t = σ(t, S t )S t dW t under a risk neutral probability measure P1 . In that case, h(t, S t ) is the value at time t of the unique self-financing portfolio process V composed of a fraction of the risk-less bond B = 1 and the risky asset S such that it replicates the European option payoff g(S T ), i.e. V T = g(S T ).

The quantity σ(t, S t ) is then interpreted as a local volatility coefficient.

In the following, we consider the uniform grid on [0, T ] given by t n i = (T /n)i, i = 0, • • • , n where n ≥ 1. We then define the following functions:

k n+ t (x) := 1 + σ(t, x) T n , k n- t (x) := 1 -σ(t, x) T n , (2.5)
λ n t (x) := 1 -k n- t (x) k n+ t (x) -k n- t (x) = 1 2 , µ n t (x) := 1 -λ n t (x) = 1 2 .
Let h n be the piecewise constant function defined on [0, T ] by

h n (t, x) = h n (t i-1 , x) if t ∈ [t n i-1 , t n i [, i ≤ n,
where the sequence (h n (t i , x)) i≤n is defined recursively by the following functional scheme:

h n (t i-1 , x) = λ n t i-1 (x)h n t i , k n- t (x)x + µ n t i-1 (x)h n t i , k n+ t (x)x , h n (T, x) = g(x), i ≤ n.
(2.6)

Our main result states that the sequence of functions (h n ) n∈N uniformly converges to h as n → ∞ under some mild conditions (called Condition D below) satisfied by the successive derivatives of h, solution to (2.4). Of course, the functional scheme (2.6) may also be considered for other choices of functions k n-and k n+ . In this paper, we restrict ourselves to the functions defined by (2.5), i.e. when the coefficients defining the binomial model are symmetric, but it is an interesting problem to study such a scheme in its full generality, i.e. for more general coefficients k n-and k n+ .

Recall that h(t, x) = Eg(S x,t (T )) where, for t ≤ T , S x,t is the unique solution to the stochastic differential equation (6.14) on the interval [t, T ] with initial condition S x,t (t) = x. We deduce that h(t, 0) = g(0). In the following, we suppose that σ is locally Lipschitz and bounded so that existence and uniqueness holds for (6.14), see for instance [START_REF] Friedman | Stochastic differential equations and applications[END_REF]Theorem 2.2,p104]. The conditions below hold in particular when the payoff function g is continuous with linear growth at infinity, e.g. g(x) = (x -K) + , K ∈ R, and the diffusion coefficient σ is smooth enough and positive as proven in [START_REF] Lépinette | Approximate hedging in a local volatility model with proportional transaction costs[END_REF]Theorem 4.1] for a more general case2 , see also [START_REF] Denis | Marchés avec coûts de transaction: approximation de Leland et arbitrage[END_REF]Section 5]. Observe that we do not suppose that the successive derivates are bounded contrarily to [START_REF] Milstein | The Probability Approach to Numerical Solution of Nonlinear Parabolic Equations[END_REF].

Condition D:

There exists a constant C > 0 independent of n such that:

|h x (t, x)| ≤ C, |h xx (t, x)| ≤ C |x| √ T -t , |h xxx (t, x)| ≤ C x 2 (T -t) , |h xt (t, x)| ≤ C T -t , |h xxt (t, x)| ≤ C |x|(T -t) 3/2 , |h tt (t, x)| ≤ C | x | (T -t) 3 2 , |h ttt (t, x)| ≤ C | x | (T -t) 5 2
,

|h xtt (t, x)| ≤ C (T -t) 2 ,
Theorem 2.1. Suppose that σ is a bounded Lipschitz function and assume that the solution to Equation (2.4) satisfies Condition (D). Then, there exists a constant C(x) > 0 depending on x, g and σ such that

sup t∈[0,T ] |h n (t, x) -h(t, x)| ≤ C|x| √ n , n > 0.
The proof is deduced from the lemmas and the corollary of Section 5. Note that, by choosing specific discretization dates, see [START_REF] Darses | Limit theorem for a modified Leland hedging strategy under constant transaction costs rate[END_REF], we should improve the convergence rate of the functional scheme. This is left for future research.

Numerical examples

In the following, we compare the performance of the functional scheme to the Crank-Nicolson finite difference method one on some examples.

The case where the diffusion coefficient only depends on time.

Suppose that σ(t, x) = σ(t), t ∈ [0, T ], only depends on time. In that case, we may show (see [START_REF] Kabanov | Mean square error for the Leland-Lott hedging strategy: convex pay-off[END_REF]) that

h(t, x) = ∞ -∞ g xe ρty-ρ 2 t 2 1 √ 2π e -y 2 2 dy, where ρ t ≥ 0 is defined by ρ 2 t = T t σ 2 u du.
Moreover, for some particular functions g, e.g. g(x) = (x -K) + , it is trivial to derive an explicit formula of h(t, x). As k n-and k n+ do not depend on the space variable, we also deduce an explicit expression of h n , solution to (2.6):

h(t n n-i , x) = 1 2 i z∈E n-i g(xz), i = 0, • • • , n, (3.7)
where the sets (E i ) i=0,••• ,n does not depend on x and are recursively defined as

E n = {1} and E n-i-1 = k n- t n n-i-1 E n-i ∪ k n+ t n n-i-1 E n-i for all i = 0, • • • , n-1.
Regarding the implementation, we may see each E i as a vector with 2 n-i components. Once these vectors computed for all, we may simulate the scheme (2.6) for any x and any arbitrary terminal condition g. Indeed, just compute the vectors G i whose components are (g(xz)) z∈E i and deduce h(t n n-i , x) as the scalar product of

G i with 1 = (1, • • • , 1) ∈ R 2 n-i .
Let us now consider the case σ(t) = σ(2 + cos(t)). Then,

ρ 2 t = 4σ 2 [(T + sin T ) -(t + sin t)] + σ 2 2 (T + 1 2 sin(2T ) -(t + 1 2 sin 2t)
and for g(x) = (x -K) + , h(t, x) is given explicitily by

h(t, x) = ∞ -∞ g xe ρty-ρ 2 t 2 1 √ 2π e -y 2 2 dy = ∞ -∞ xe ρty-ρ 2 t 2 -K + 1 √ 2π e -y 2 2 dy = ∞ ρ -1 t ln( K x )+ ρ 2 t 2 xe ρty-ρ 2 t 2 -K 1 √ 2π e -y 2 2 dy = ∞ ρ -1 t ln( K x )+ ρ 2 t 2 x e -(y-ρ t ) 2 2 √ 2π dy -K 1 -Φ ln( K x ) + ρ 2 t 2 ρ t = x 1 -Φ ln( K x ) + ρ 2 t 2 ρ t -ρ t -KΦ ln( x K ) - ρ 2 t 2 ρ t = xΦ ln( x K ) + ρ 2 t 2 ρ t -KΦ ln( x K ) - ρ 2 t 2 ρ t .
Therefore, we are in position to numerically compare the convergence error between the explicit solution given above and the functional scheme as well as the finite difference approximations. We first consider the parameters σ = 0.5, T = 10, K = 100, x = 120 and 30 discrete dates. As we may observe in Figure 1, with a space discretization of 400 points, the finite difference method does not provide a good approximation and, moreover, the monotonicity we expect to is not satisfied. With the parameters σ = 1.5, T = 1, K = 115, x = 95 and the same discretization in time and in space (see Figure 22), the finite difference method provides an approximating function which explodes as time tends to 0.

At last, with the parameters σ = 0.1, T = 5, K = 115, we consider 25 discrete dates and 400 points for the discretization in time. We observe that the approximating function given by the finite difference method is not very efficient as well.

In order to better compare the two methods, we have computed the needed time for several strikes to obtain a reasonable convergence error for the value of h(0, x) with respect to the analytic expression of the later. Initially, the time discretization is done with 10 dates while the space is discretized in 400 points. If the error is larger than 5%, we add as many discrete dates as necessary (see Tables 4, 5, 6 and6). Clearly, the results are bad for the finite difference method, see Table 5. With the parameters σ = 1.5, T = 1, x = 95, we get the following: Notice that the results given by the finite difference method are not ac-ceptable, e.g. when K is too large: 3.2. The case where the diffusion coefficient only depends on space.

Let us consider the diffusion coefficient σ(t, x) = σ(x) = σ * (1 + e -x 2 ). As h(t, x) does not admit any analytic expression, we first evaluate it by using its probabilistic representation through a Monte-Carlo discretization of the associated diffusion process (6.14). The discretization in time is composed of 10000 dates with 50000 samples for the Monte-Carlo method and 30 dates for the two other methods while we discretize the space variable in 400 points.

As we may observe in Figure 8, the finite difference method does not provide a non increasing function in time as expected and the approximation is not accurate. In Table 11, we see that refining more the space variable grid does not significantly improve the results while the computation time increases. With the parameters σ = 1.5, T = 1, K = 115, x = 95 and only 25 discrete dates, the function scheme provides a rather good approximation contrarily to the finite difference method: At last, we consider the parameters σ = 2, T = 5, K = 10, x = 95 with only 20 discrete dates for the functional scheme. We get satisfactory results in less that 0.26 seconds while with 30 discrete dates for the finite difference methods, the computation time is large. In particular, even if the finite difference method provides acceptable results with only 10 discrete dates with an error less than 5% it takes more than 3.60 secondes to computes them.

The case where the diffusion coefficient depends on time and space.

We consider the volatility function σ(t, x) = σ(1 + 0.01x 2 1+0.01x 2 + t). The number of discrete dates is 25 for the functional scheme as well as the finite difference method and 400 points for the space variable. The Monte Carlo discretization is based on 10000 discrete dates and 50000 trajectories. When the parameters For the parameters σ = 0.1, T = 5, K = 110, x = 95, the functional scheme over evaluates the Monte Carlo method while monotonicity is not satisfied for the finite difference method. Notice that it is possible to execute parallel computations by means of programming interfaces (e.g. OpenMP or MPI) i.e. we may use several calculus units simultaneously and reduce computation time. Using n ≥ 1 processors, we should theoretically divide the needed time by n. In practice, that depends on how the executed parallel calculus are interconnected. In the finite difference method, when computing the value function at some point of the grid, we need to estimate some expressions depending on the adjacent nodes that may be calculated and used by other processors. This should increase the computation time, as it is necessary to wait for some computations to be executed before starting new ones. On the contrary, for the functional scheme, we do not face this problem as the nodes of the binomial tree at the same level are independent. This is confirmed by the realized time savings we observe to be very closed to the theoretical ones (e.g. 50% with two cores and 75% with 4 cores), see Table 21 with the example of the current subsection and the parameters σ = 0.5, T = 10, K = 100, S 0 = 120. Therefore, even with parallel computing, the functional scheme should outperform the finite difference method in the examples we have presented. Let us consider the payoff function h(x) = (K -x) + . The case when h(x) = (x-K) + may be deduced by the Call-Put parity (x-K) + -(K -x) + = x-K.

When the volatility function only depends on time (see Section 3.1), we obtain very good approximations, see Figure 22 and Table 23, when the space variable interval is [0, 10 4 ] discretised in 2000 points. This is confirmed by other sets of parameters.

The case of an unbounded payoff function.

Let us consider the payoff function h(x) = (x-K) 8 and let us fix the following parameters: σ = 0.25, T = 1 and K = 0.5. We then compute t → h(t, x) for x = 1. We compare the obtained results to the approach based on the probabilistic representation of h(t, x), i.e. by implementing the Monte Carlo (MC) technique with a time discretization of 50000 dates and a sample of 100000 trajectories. First, we observe that the (FE) method is very performant when σ(t, x) = σ is a constant function and the interval space [0, 10] is discretised in 500 points. To see it, we estimate the relative error of the (FE) approximation with respect to the (MC) method, see In this example, the functional scheme seems to be more unstable even if it only takes 3 m.s. to get an error less than 5%. Nevertheless it takes too much time to obtain an error less than 1%, see Table 4.2.

We now consider the volatility function

σ(t, x) = σ 1 + x 2 100 1 + x 2 100 + t .
In that case, it is difficult to obtain the desired convergence as the result is very sensitive to the chosen space variable and its discretization as well. It appears that we reduce the error by increasing the discretization in time but it is not possible to lower the error as much as desired, see Table 26.

Surprisingly, if we reduce the number of discrete points of the space variable to 100 points, the error is smaller as observed in Table 27. Nevertheless, even by increasing the number of discrete dates to 10 7 and more, it is not possible to reduce as much as desired the convergence error. If we enlarge the space variable to [0, 25] the error is 837% for 2000 discrete dates so that we need to increase more the number of discrete dates. On the contrary, the functional scheme appears to be efficient, see Table 28. In particular, it is possible to reduce as much as desired the convergence error. It takes about 145 seconds to obtain an error less that 5%. Moreover, a code parallelisation should significantly reduce the needed computation time.

Proof

From a line to the next one, we shall use the same notation C for some distinct constants independent of n. Also, notice that there is a constant

M > 0 such that |k n- t | + |k n+ t | ≤ M for all n and t ∈ [0, T ]. We also use the notation ∆t n i = t n i -t n-1 i = T /n, i ≤ n.
At last, observe that the results below are trivial when x = 0 hence we only provide the proofs when x = 0. Lemma 5.1. Suppose that the conditions of Theorem (2.1) hold. Then, there exists a constant C > 0 which does not depend on x such that

sup t∈[t n n-1 ,T ] |h n (t, x) -h(t, x)| ≤ C|x| √ n .
Proof. Let us define δ n t =| h(t, x)-h n (t, x) |. Recall that h(t, x) = Eg (S x,t (T )) where S x,t is the solution to (6.14). We have:

δ n t = | h(t, x) - 1 2 g(k n- t n n-1 x) - 1 2 g(k n+ t n n-1 x) | ≤| h(t, x) -g(x) | + 1 2 | g(k n- t n n-1 x) -g(x) | + 1 2 | g(k n+ t n n-1 x) -g(x) | ≤ sup t n n-1 ≤t≤T | h(t, x) -g(x) | + L g 2 (1 -k d t n n-1 ) | x | + L g 2 (k u t n n-1 -1) | x | ≤ L g sup t n n-1 ≤t≤T E | S x,t (T ) -x | + L g 2 σ(t, x) T n | x | + L g 2 σ(t, x) T n | x | ≤ C|x| T n + L g 2 σ * T n | x |,
where the last inequality is deduced from Lemma 6.1 and

σ * = sup t∈[0,T ],x∈R |σ(t, x)|. 2
Lemma 5.2. Suppose that the conditions of Theorem (2.1) hold. Then, there exists a constant C > 0 such that

|h n (t n n-2 , x) -h(t n n-2 , x)| ≤ C|x| √ n .
Proof. Let us introduce δ n t-2 =| h n (t n n-2 , x) -h(t n n-2 , x) |. We have

δ n t-2 = | 1 2 h n (t n n-1 , k n- t n n-2 x) + 1 2 h n (t n n-1 , k n+ t n n-2 x) -h(t n n-2 , x) | ≤ 1 2 | h n (t n n-1 , k n- t n n-2 x) -h(t n n-2 , x) | + 1 2 | h n (t n n-1 , k n+ t n n-2 x) -h(t n n-2 , x) | ≤ 1 2 | h n (t n n-1 , k n- t n n-2 x) -h(t n n-1 , k n- t n n-2 x) | + 1 2 | h(t n n-1 , k n- t n n-2 x) -h(t n n-2 , k n- t n n-2 x) | + 1 2 | h(t n n-2 , k n- t n n-2 x) -h(t n n-2 , x) |) + 1 2 (| h n (t n n-1 , k n+ t n n-2 x) -h(t n n-1 , k n+ t n n-2 x) | + 1 2 | h(t n n-1 , k n+ t n n-2 x) -h(t n n-2 , k n+ t n n-2 x) | + 1 2 | h(t n n-2 , k n+ t n n-2 x) -h(t n n-2 , x) |).
Under Condition D, by the mean value theorem, we deduce that

| h(t n n-2 , x) -h(t n n-2 , k n- t n n-2 x) | ≤ C | x -k n- t n n-2 x |≤ C | x | . | 1 -k n- t n n-2 | ≤ C | x | T n , | h(t n n-2 , x) -h(t n n-2 , k n+ t n n-2 x) | ≤ C | x -k n+ t n n-2 x |≤ C | x | . | 1 -k n+ t n n-2 | ≤ C | x | T n .
Note that, under Condition D, we have

| h t (t, x) |= x 2 σ 2 (t, x) 2 | h xx (t, x) |≤ C | x | √ T -t . (5.8)
Therefore, by the mean value theorem, since k n-

t n n-2 , k n+ t n n-2 ≤ M , | h(t n n-1 , k n- t n n-2 x) -h(t n n-2 , k n- t n n-2 x) | ≤ CM | x | T -t n n-1 ∆t n n-1 ≤ CM | x | T n , | h(t n n-1 , k n+ t n n-2 x) -h(t n n-2 , k n+ t n n-2 x) | ≤ CM | x | T -t n n-1 ∆t n n-1 ≤ CM | x | T n .
At last, by Lemma 5.1,

γ n -(t n n-1 , x) :=| h n (t n n-1 , k n- t n n-2 x) -h(t n n-1 , k n- t n n-2 x) | ≤ Ck n- t n n-2 |x| T n , γ n + (t n n-1 , x) :=| h n (t n n-1 , k n+ t n n-2 x) -h(t n n-1 , k n+ t n n-2 x) | ≤ Ck n+ t n n-2 |x| T n . As k n- t n n-2 + k n+ t n n-2 = 2, we deduce that 1 2 γ n -(t n n-1 , x) + 1 2 γ n + (t n n-1 , x) ≤ C|x| T n .
The conclusion follows.2

Lemma 5.3. Suppose that the conditions of Theorem (2.1) hold. Then, there exists a constant C > 0 such that

|h n (t n i , x) -h(t n i , x)| ≤ C|x| √ n , for all i ≤ n -2.
Proof. We first prove the result for i = n -2 and, then, we generalize the result by induction. Let us introduce the function

F i (x) = 1 2 h(t n i , k n- t n i-1 x) + 1 2 h(t n i , k n+ t n i-1 x)
. By the Taylor formula, we get that

F i (x) = 1 2 h(t n i-1 , x) + R 11 i (x) + R 12 i (x) + R 13 i (x) + 1 2 h(t n i-1 , x) + R 21 i (x) + R 22 i (x) + R 23 i (x) .
where, for some constants

α i , i = 1, • • • , 4, and tn i-1 ∈ [t n i-1 , t n i ], x = x + θx(k n- t n i-1 -1), θ ∈ [0, 1], tn i-1 ∈ [t n i-1 , t n i ], x = x + θx(k n+ t n i-1 -1), θ ∈ [0, 1], we have R 11 i (x) = h t (t n i-1 , x) T n -h x (t n i-1 , x)xσ(t n i-1 , x) T n , R 12 i (x) = -h tx (t n i-1 , x)xσ(t n i-1 , x) T n 3 2 + 1 2 h xx (t n i-1 , x)x 2 σ 2 (t n i-1 , x) T n + 1 2 h tt (t n i-1 , x) T n 2 , R 13 i (x) = α 1 h ttt ( tn i-1 , x) T n 3 -α 2 h ttx ( tn i-1 , x)σ(t n i-1 , x)x T n 5/2 +α 3 h txx ( tn i-1 , x)σ 2 (t n i-1 , x)x 2 T n 2 -α 4 h xxx ( tn i-1 , x)σ 3 (t n i-1 , x)x 3 T n 3/2 , R 21 i (x) = h t (t n i-1 , x) T n + h x (t n i-1 , x)xσ(t n i-1 , x) T n , R 22 i (x) = h tx (t n i-1 , x)xσ(t n i-1 , x) T n 3 2 + 1 2 h xx (t n i-1 , x)x 2 σ 2 (t n i-1 , x) T n + 1 2 h tt (t n i-1 , x) T n 2 , R 13 i (x) = α 1 h ttt ( tn i-1 , x) T n 3 + α 2 h ttx ( tn i-1 , x)σ(t n i-1 , x)x T n 5/2 +α 3 h txx ( tn i-1 , x)σ 2 (t n i-1 , x)x 2 T n 2 + α 4 h xxx ( tn i-1 , x)σ 3 (t n i-1 , x)x 3 T n 3/2 . We deduce that F i (x) = h(t n i-1 , x) + R 1 i (x) + R 2 i (x) + R 3 i (x) where R 1 i (x) = 1 2 (R 11 i (x) + R 21 i (x)), R 2 i (x) = 1 2 (R 12 i (x) + R 22 i (x)) and R 3 i (x) = 1 2 (R 13 i (x) + R 23 i (x)
). As h is the solution to (2.4), we may simplifly F i (x) as

F i (x) = h(t n i-1 , x) + 1 2 h tt (t n i-1 , x) T n 2 + R 3 i (x).
Under Condition D, there exists a constant C independant of i and n such that

1 2 h tt (t n i-1 , x) ≤ C | x | (T -t n i-1 ) 3 2
.

(5.9)

On the other hand, since tn i-1 , tn i-1 ∈ [t n i-1 , t n i ] and |x| + |x| ≤ C|x| where C > 0 is independent of i and n, the residual error R 3 i (x) is the sum of terms which may be dominated under Condition D using the following inequalities:

α 1 h ttt ( tn i-1 , x) ≤ C | x | (T -t n i ) 5 2
.

(5.10)

α 2 h xtt ( tn i-1 , x)σ(t n i-1 , x)x ≤ C | x | (T -t n i ) 2 .
(5.11)

α 3 h txx ( tn i-1 , x)σ 2 (t n i-1 , x)x 2 ≤ C | x | (T -t n i ) 3 2
.

(5.12)

The inequalities above are also satisfied if we replace tn i-1 by tn i-1 and x by x. Recall that the constant C may denote distinct constants that change for a line to the next one but this constant does not depend on i and n. At last, by the Taylor formula, we rewrite the last term as

i = α 4 h xxx ( tn i-1 , x)σ 3 (t n i-1 , x) -h xxx ( tn i-1 , x)σ 3 (t n i-1 , x) x 3 T n 3/2 , = α 4 h xxxt ( tn i-1 , x)x 3 ( tn i-1 -tn i-1 ) T n 3/2 + α 4 h xxxx ( tn i-1 , x)x 3 (x -x) T n 3/2
, where tn

i-1 ∈ [ tn i-1 , tn i-1 ] and x = αx + (1 -α)x for some α ∈ [0, 1]. It is easily seen that | tn i-1 -tn i-1 | ≤ CT /n
for some constant C independent of i and n. Similarly, by definition of x and x, there exists a constant C such that |x -x| ≤ C|x| T /n. Therefore,

| i | ≤ C | x | (T -t n i ) 2 T n 5/2 + C | x | (T -t n i ) 3/2 T n 2 . (5.13) Let us introduce δ n t i-1 (x) := h n (t n i-1 , x) -h(t n i-1 , x) , i ≤ n -2.
From above, we get that

δ n t i-1 (x) = 1 2 h n (t n i , k n- t n i-1 x) + 1 2 h n (t n i , k n+ t n i-1 x) -h(t n i-1 , x) = 1 2 h n (t n i , k n- t n i-1 x) + 1 2 h n (t n i , k n+ t n i-1 x) -F i (x) + 1 2 h tt (t n i-1 , x) T n 2 + R 3 i (x) ≤ 1 2 δ n t i (k n- t n i-1 x) + 1 2 δ n t i (k n+ t n i-1 x) + 1 2 T n 2 h tt (t n i-1 , x) + R 3 i (x) .
Using the inequalities (5.9), (5.10),• • • , (5.13), we then deduce a constant C > 0 independent of i and n such that

δ n t n i-1 (x) ≤ 1 2 δ n t i (k n- t n i-1 x) + 1 2 δ n t i (k n+ t n i-1 x) + C|x| (T -t n i ) 1 2
T n

3 2 + C|x| (T -t n i ) 3 2 T n 2 + C|x| (T -t n i ) 2 T n 5 2 + C|x| (T -t n i ) 5 2
T n

3

.

Recall that, by Lemma 5. T n Both with Lemma 5.1 and 5.2 and the inequalities above, we may conclude. We then conclude by Lemma 5.3.2

Appendix

Let us consider the unique solution S x,t , t ∈ [0, T ], to the stochastic differential equation dS t,x (u) = S t,x (u)σ(u, S t,x (u))dW u , u ∈ [t, T ], S t,x (t) = x ∈ R, (6.14) where W is a standard Brownian motion. Applying the Gronwall lemma, we deduce that:

E | S x,t (r) | 2 ≤ 8x 2 exp (8σ * (T -t)) ≤ Cx 2 ,
where C does not depend on x. By the Ito isometry, we then deduce that E(S x,t (T ) -x) 2 = E T t σ 2 (u, S x,t (u))S 2 x,t (u)du Using the inequality above, we may conclude.2

Fig 1 .

 1 Fig 1. Numerical example with σ = 0.5, T = 10, K = 100, x = 120.

Fig 4 .

 4 Fig 4. Values at t = 0 by the functional scheme with parameters σ = 0.5, T = 10, x = 120.

Fig 5 .

 5 Fig 5. Values at t = 0 by finite differences with parameters σ = 0.5, T = 10, x = 120.

Fig 2 .

 2 Fig 2. Numerical example with σ = 1.5, T = 1, K = 115, x = 95.

Fig 3 .

 3 Fig 3. Numerical example with σ = 0.1, T = 5, K = 115, x = 115.

Fig 6 .

 6 Fig 6. Values at t = 0 by the functional scheme with parameters σ = 1.5, T = 1, x = 95.

Fig 7 .

 7 Fig 7. Values at t = 0 by the functional scheme with parameters σ = 1.5, T = 1, x = 95.

Fig 8 .

 8 Fig 8. Numerical example with σ = 0.5, T = 10, K = 100, x = 120.

Fig 10 .

 10 Fig 10. Values at t = 0 by the finite differences with parameters σ = 0.5, T = 10, x = 120.

Fig 11 .

 11 Fig 11. Values at t = 0 by the finite differences with parameters σ = 0.5, T = 10, x = 120 and K = 100 when refining the space variable grid.

Fig 12 .

 12 Fig 12. Numerical example with σ = 1.5, T = 1, K = 115, x = 95.

Fig 13 .

 13 Fig 13. Values at t = 0 by the functional scheme with parameters σ = 1.5, T = 1, x = 95.

Fig 14 .

 14 Fig 14. Numerical example with σ = 2, T = 5, K = 10, x = 95

Fig 15 .

 15 Fig 15. Numerical example with σ = 0.5, T = 10, K = 100, x = 120.

Fig 17 .

 17 Fig 17. Values at t = 0 by finite differences with parameters σ = 0.5, T = 1, x = 120.

Fig 18 .

 18 Fig 18. Numerical example with σ = 0.1, T = 5, K = 110, x = 95.

Fig 20 .

 20 Fig 20. Values at t = 0 for finite differences with parameters σ = 0.1, T = 5, x = 95.

Fig 21 . 4 . 4 . 1 .

 21441 Fig 21. Time savings with the parallel functional scheme when σ = 0.5, T = 10, K = 100, S 0 = 120.

Fig 22 .

 22 Fig 22. Numerical example with σ = 1.5, T = 1, K = 115, x = 95.

Fig 23 .

 23 Fig 23. Values at t = 0 by the FEM with parameters σ = 1.5, T = 1, K = 115, x = 95.

Fig 28 .

 28 Fig 28. Approximation of h(0, x) with the (FS) method and non constant volatility.

3 .

 3 Since (T -u) -1 ≥ (T -t n j ) -1 if u ∈ [t n j , t n j+1 ],we deduce that

  Values at t = 0 by the functional scheme with parameters σ = 0.5, T = 1, x = 120.

	Strike Result Error Time Discretization Computation time (s) 50 83.21 4.36% 11 < 0.01 60 78.26 4.52% 11 < 0.01 70 73.64 3.66% 11 < 0.01 80 69.61 4.24% 11 < 0.01 90 66.18 3.41% 11 < 0.01 100 63.15 4.77% 11 < 0.01 110 60.32 4.04% 11 < 0.01 Error Time Discretization Computation time (s) 50 0.02 -99.98% 30 11.42 60 0.01 -99.99% 30 13.00 70 0.00 -100.00% 30 11.82 80 0.00 -100.00% 30 12.23 90 0.00 -100.00% 30 11.77 100 0.00 -100.00% 30 12.42 Fig 16. Strike Result 110 0.00 -100.00% 30 11.77

Table 24 .

 24

	number of dates	(FE) method	Time (s)	Error
	5	1.1201	0.001	215.32%
	10	0.5969	0.001	68.02%
	20	0.4508	0.002	26.91%
	50	0.3855	0.005	8.52%
	100	0.3666	0.009	3.19%
	200	0.3576	0.017	0.66%
	300	0.3547	0.025	0.16%

Fig 24. Approximation of h(0, x) with the (FE) method.

  Approximation of h(0, x) with the (FS) method and constant volatility.

	Number of dates	Functional Scheme	Time (s)	Error
	5		0.1831	0.000	93.96%
	6		0.3331	0.000	6.66%
	7		0.3344	0.000	6.23%
	8		0.2280	0.000	55.81%
	9		0.2382	0.000	49.15%
	10		0.3378	0.001	5.16%
	11		0.2543	0.001	39.70%
	12		0.2608	0.001	36.23%
	13		0.3402	0.003	4.41%
	14		0.3409	0.006	4.21%
	15		0.3415	0.013	4.03%
	16		0.2799	0.012	26.91%
	17		0.2835	0.024	25.29%
	18		0.2868	0.049	23.87%
	19		0.3433	0.197	3.48%
	20		0.2925	0.196	21.47%
	21		0.2949	0.400	20.44%
	22		0.3443	1.571	3.17%
	23		0.3446	3.135	3.08%
	24		0.3449	6.285	3.00%
	25		0.3031	6.237	17.20%
	26		0.3454	24.917	2.86%
	27		0.3456	51.146	2.80%
	28		0.3458	107.718	2.74%
	29		0.3460	213.635	2.68%
	Fig 25. nbIterations	Finite Elements Method	Time (s)	Error
	5		110,8036	0.001	644.37%
	10		63.7537	0.001	328.29%
	20		39.5068	0.003	165.40%
	50		26.7455	0.006	79.67%
	100		23.0440	0.011	54.81%
	200		21.3239	0.022	43.25%
	300		20.7713	0.034	39.54%

Fig 26. Approximation of h(0, x) with the (FE) method and non constant volatility

2

  Corollary 5.4. Suppose that the conditions of Theorem (2.1) hold. Then, there exists a constant C > 0 such thatProof. Let t ∈ [0, T ] be such that t ∈ [t n i-1 , t n i )for some i ≥ 1. By Lemma 5.1, we may suppose that i ≤ n -2. Then,|h n (t, x) -h(t, x)| = |h n (t n i-1 , x) -h(t, x)| ≤ |h n (t n i-1 , x) -h(t n i-1 , x)| + |h(t n i-1 , x) -h(t, x)|.By the mean value theorem and Inequality (5.8), since T -t n i ≥ T /n, we deduce that

	sup t∈[0,T ]	|h n (t, x) -h(t, x)| ≤	C|x| √ n	.
	|h(t n i-1 , x) -h(t, x)| ≤	C|x| T -t n i	∆t n i ≤	C|x| √ n	.

  Lemma 6.1. Suppose that t ∈ [0, T ]. Let S x,t be the solution to the stochastic differential equation(6.14). If σ is bounded by a constant σ * > 0, there existsa constant C such that (u) ≤ Cx 2 , E(S t,x (T ) -x) 2 ≤ Cx 2 (T -t).Proof. By the Doob's inequality, we obtain that for every t ≤ r ≤ T :φ(r) := E | sup t≤u≤r S x,t (u) | 2 ≤ 4E | S x,t (r) | 2As S x,t (r) = x + r t σ(u, S x,t (u))S x,t (u)dW u , using the inequality (a + b) 2 ≤ 2(a 2 + b 2 ) and the Ito isometry, we get that :

			r	2
	E sup t≤u≤T t,x φ(r) ≤ 8x 2 + 8E S 2 ≤ 8x 2 + 8E ≤ 8x 2 + 8(σ * ) 2 ≤ 8x 2 + 8(σ * ) 2	t t	σ(u, S x,t (u))S x,t (u)dW u σ 2 (u, S x,t (u))S 2 x,t (u)du r t r ES 2 x,t (u)du r t φ(u)du.

If we suppose that the risk-free interest rate is r = 0

Here, the coefficient γ n of[START_REF] Lépinette | Approximate hedging in a local volatility model with proportional transaction costs[END_REF] appears to be a constant independent of n.


Héctor Alarcón

Jean-Christophe Géminard

Francisco Melo

Patterning of a cohesionless granular layer under pure shear

Keywords: 89.75.Kd: Pattern formation in complex systems, 83.60.Uv: Rheology: fracture, 45.70.Qj: Pattern formation in granular matter

The response of a thin layer of granular material to an external pure shear imposed at its base is investigated. The experiments show that, even for non-cohesive materials, the resulting deformation of the material is inhomogeneous. Indeed, a novel smooth pattern, consisting in a periodic modulation of the shear deformation of the free surface, is revealed by an image-correlation technique. These observations are in contrast with the observation of the fracture pattern previously observed in cohesive granular materials subjected to stretching. For cohesive materials, the instability is due to the weakening of the material which results from the rupture of capillary bridges that bond the grains to one another. For non-cohesive materials, the rupture of the capillary bridges cannot be invoked anymore. We show that the instability results from the decrease of friction upon shearing.

I. INTRODUCTION

Cohesive granular materials are characterized by a network of liquid bridges responsible of attractive capillary forces between particles [1][2][3]. Various regimes of cohesion are identified depending on the liquid content, leading to different scalings for the cohesion force [2,4,5]. Regardless of fluid content, a common feature of cohesive material is a weakening due to a decrease of both the associated adhesion force when a single bridge elongates [6] and the overall number of bonds which collapse when excessively stretched [7]. This effect, which we referred to as strain softening [8], is observed above a critical stress associated with the cohesion due to the capillary bridges at grain contacts [9] and is responsible for the relatively low range of plasticity in cohesive granular materials.

In a recent work [8] we explored the tensile response of a horizontal layer of cohesive material subjected to homogeneous deformation in its bottom plane. We showed that "strain softening" is responsible for the nearly periodic modulation of the strain field along the pulling axis that develops as soon as the external deformation is imposed. The associated wavelength depends linearly on the layer thickness, is almost independent of particle size and depends linearly on the relative humidity. The flexural deformation of a cohesive granular layer reveals similar features [10]. To establish a more fundamental connection between the pattern features and the intrinsic properties of the granular material, we developed experimental methods for the assessment of the cohesion and shear modulus as function of the particles size and relative humidity [11].

Here, we explore the response of a granular layer to pure shear at its base. For important cohesion, we observe that cracks appear and form a periodic network. However, when cohesion is reduced, the layer deforms without fracturing and a pattern, consisting in rather a smooth, periodic, modulation of the layer thickness is observed, instead. For vanishing cohesion, the pattern does not disappear but the amplitude of modulation is so small that it is only revealed by the use of image correlation techniques. A simple model including, besides cohesion, the decrease of friction due to the induced dilation of the material accounts for the typical size of the structure.

II. EXPERIMENTAL SETUP AND PROTOCOLS

The experiment consists in imposing an in-plane deformation at the base of a thin layer of granular material. The mechanical system is composed of four linear actuators (Thorlabs Z825BV), placed in the xy axes (Fig. 1) (see Ref. [11] for details). Actuators are attached to four plexiglass blocks respectively, each of them holding one arm of a latex membrane cut in cross-shape (thickness 0.5 mm, width 40 cm). The membrane leans on a horizontal table that prevents bending due to weight and ensures planar deformation. Pure shear deformation at center is obtained by stretching two opposite arms and by shortening the arms in the perpendicular direction at the same velocity. We checked the amplitude and homogeneity of the deformation through digital image correlation.

The sample is prepared by pouring granular material into a circular mold of given height h (from 1 to 10 mm, to within 0.1 mm) and internal diameter 7 cm, leaning on the membrane, at center. The excess of grains is gently removed through the horizontal displacement a rod. After removal of the mold, we obtain a disk of granular material of well-defined thickness h. We used either glass or brass spherical-beads, the large density of brass making possible to explore the vanishing small cohesion regime, the weight of the grains overcoming more easily the capillary forces. Two ranges of particle diameter d are considered for both materials: (0 -45) or (150 -200) µm for glass, and (75 -106) or (212 -300) µm for brass.

The whole experimental device is placed in a chamber whose atmosphere is maintained at constant relative humidity, R H [equilibrated with saturated salt solutions and monitored using a humidity meter (Lutron HT-3015)]. Prior to any experimental run, particles are sonicated in acetone, rinsed in pure water and dried, to prevent undesired cohesion due to contamination. Cohesion is measured, in addition, as described in Ref. [11].

The free surface of the sample is imaged from above by using a digital camera (Nikon DMX1200). An annular light-source (ring of LEDs, Fig. 1) placed 1 cm above the sample, at center, provides the grazing lighting adequate to observation of tiny displacements at the layer surface. The typical experimental run consists in applying the pure-shear deformation at the membrane in quasi-static manner and in determining precisely the displacement fields at the free surface of the granular layer. In practice, the imposed shear strain, θ is increased by constant steps ∆θ, of the order of 10 -3 . After each step, the free surface at rest is imaged with a resolution of (3840 × 3072) px 2 . In order to extract the resulting displacement field, the following correlation procedure is then applied to the successive images: We define a sliding window (typical surface area 1.5 mm 2 , containing enough particles to be considered as a coarse grain), and scan the whole image (typical displacement 0.2 mm). The procedure achieves a spatial resolution of the order of 1 mm and a resolution of the order of a few µm in the local grain displacement. Shear and vertical vorticity fields are obtained through differentiation of displacement field [11].

Grazing light Motors

III. EXPERIMENTAL RESULTS

When the thickness, h, and the cohesion, σ s , are sufficiently large, we observe fracturing of the granular layer (Fig. 2). The fractures organize in rather periodic pattern. The typical distance between neighbour fractures, λ, increases with both the layer thickness and cohesion. It is of particular interest that a pattern is still present even when fracturing is not observed in Fig. 2 (h = 3 mm and R H = 39 %, for instance). This first observation is confirmed by the results obtained with larger grains (Fig. 3): even if not visible on simple images of the free surface, the squared patterns is revealed by the shear or vorticity fields, whatever the thickness and relative humidity in the experimental range (Fig. 4). In addition, note that, as expected, the typical distance between neighboring fractures is no longer a function of R H in the limit of vanishing cohesion (Fig. 3b). The previous results reveal that, upon shear, the granular layer is subjected to an instability leading to a modulation of the in-plane strain field. A sensitive method to determine the onset of instability consists in considering the amplitude of the modulation of the grain displacement at the free surface as function of the shear-strain imposed in the bottom plane. In order to measure a representative averaged amplitude of the modulation, we select a line profile along a perpendicular to the fractures. The displacement along the fracture, ∆d, is then measured revealing a well-defined modulation at the pattern wavelength, λ. In Fig. 5, we report the averaged amplitude of such modulation, < ∆d >, as a function of θ. We observe that, < ∆d > is nearly zero for small θ, but significantly increases for θ above a critical value. In addition, we note that the maximum of the averaged amplitude < ∆d > increases the layer thickness, h. Even if the behavior of the onset as function of h is difficult to assess from our experimental data sets, the onset of instability is clearly smaller for larger cohesion, at a given h.

✓ = 0.083 ✏ xy ! z h = 3 mm h = 7 mm
In order to explore the regime of vanishing cohesion in the same range of grain size, we now use brass particles and reduce the relative humidity as much as possible. In a dry system, the cohesion originates from the Van der Waals attraction force between the particles. We estimate from the Hamaker's constant, A H ≈ 4 × 10 -19 J, and from the typical size of the asperities at the particles surface, D ≈ 200 nm (obtained with an atomic force microscope [11]), the tensile stress associated with an individual asperity, separated from a flat surface by a distance z of atomic scale (z = z a ≈ 1nm) is about A H D/(3d 2 z 2 a ) ≈ 0.02 Pa. Considering that cohesion is due to the contribution of all asperities, of surface density ρ s , in the region of contact, assuming the separation z = z a + 2r 2 /d at a distance r from the center of the contact region, we get σ s ≈ πA H Dρ s /(3dz a ). From [11], we estimate that the distance between asperities is of about 5D, such that ρ s ≈ 1/(5D) 2 and, thus, σ s ≈ 1 Pa. This value is much smaller that the typical pressure P = ρgh ≈ 400 Pa (with ρ b = 8 × 10 3 kg/m 3 and h ≈ 0.5 cm) at the base of the granular layer. When shear is imposed to a such non-cohesive granular material, the cellular pattern is revealed by image correlation analysis, only. The vorticity fields are the most suitable to reveal the structure (Fig. 6). We again observe that the wavelength, λ, increases with the layer thickness, h. We also note that λ increases with the grain size.

IV. ANALYSIS

Our analysis of the system behavior is based on the simplest heuristic model accounting for friction and, eventually, cohesion.

When a dry granular material is deformed, the main energetic cost is from the solid friction between the grains [8, 10,11]. The simplest way to account for friction is to write that, when two surfaces in regard are displaced with respect to one another, the energetic cost is proportional to the friction coefficient µ and local pressure P . However, the friction is known to decrease upon shear due to the dilation of the material, an assumption that is well supported by early works [12]. Starting from the initial value µ s previous to shear, the friction coefficient is decreased by ∆µ for a typical displacement δ (of the order of the grain radius) of the two surfaces in regard. Then, the friction coefficient remains constant and equal to µ d ≡ µ s -∆µ. Here, seeking for simplicity, we will assume that the frictional coefficient is constant and equal to µ = µ d once the material is sufficiently deformed, but that an additional energetic cost, Γ = ∆µP δ, is to be paid to initially deform the material.

In a dry granular material, the local pressure P is due to the weight of the grains, only. We have, at a distance z from the free surface, P = ρgz. In order to account for the results obtained for small, but not negligible cohesion, we generalize slightly the model by introducing the energetic cost associated with the liquid bridges. To do so, we take into account that the liquid bridges are responsible for a tensile pressure σ s , independent of the depth z, which must be added to the pressure P . However, when the material is deformed, σ s works only on a distance d b that is of the order of the bridge size. For the humidity content considered here, the typical size of the capillary bridge is of the order of the size of the asperities at the grain surface. Thus, d b δ. When the relative displacement of two grains exceeds δ, the bridges are broken. Considering that the force associated with the tensile pressure only works over the distance d b , we obtain that the contribution of the cohesion to Γ is of the order of (µ + ∆µ)σ s d b . In the following, we thus consider that Γ = ∆µP δ + µσ s d b , where we considered ∆µ µ. Let us assume that, upon uniform shear in the bottom plane, the granular layer forms a series stripes of width λ that are not deformed, but separated from the membrane and from one another by shear bands. The energetic cost associated with the deformation is thus the sum of the cost of the shear band in the bottom plane, plus the cost of the vertical shear bands. At bottom, we have:

E h = Γ h λ + µ ρgh +λ/2
-λ/2 θx dx. The first contribution Γ h λ is the cost of the fracture at the bottom which is the sum of the excess energy until the liquid bridges are broken, µσ s d b , and of the excess energy to initiate the shear deformation, ∆µ ρgh δ. The second contribution is that of friction, the local pressure being P = ρgh. Considering the dimensionless cost per unit length E h ≡ E h /(µ ρgh λ), we write:

E h = d b + ∆µ µ δ + θ 4 λ. (1)
where we introduced ≡ σ s /(ρgh). In the same way, at the vertical walls between the bands, we have:

E v = Γ v h + µ ρg [ h 0 (h -z) dz]
θλ where the cost of the vertical wall is Γ v = µσ s d b h + ∆µ ρg(h/2) δh (the pressure P is averaged over the thickness). The dimensionless cost per unit length is:

E v = 1 2 2 d b λ + ∆µ µ δ λ + θ h. (2)
The sum E h + E v is the total cost (i.e. the energy lost by friction) per unit length in the volume upon deformation of the granular layer.

In order to determine if the layer is unstable (with the respect of the modulation with the wavelength λ) and to obtain the most unstable wavelength, we first estimate the energy E s that would be elastically loaded in the system, supposed to deform homogeneously, for the same imposed shear at the bottom. The energy of this homogeneous state is obtained, provided the knowledge of the shear modulus G of the granular layer, by writing

E s = λ h 0 1 2 G θ 2 dz.
In the absence of cohesion, G = α P where α = Ed/(6Jl R ) [11]. The prefactor α thus depends on physical properties of the grains, such as the Young modulus E, the yield stress J and the typical size of the surface asperities l R . We have:

E s = 1 4 α µ (1 + 2 ) h θ 2 , (3)
where the coefficient µ appears only because our choice of the energy scale (i.e. µ ρghλ).

In order to determine the onset of the instability and the most unstable wavelength, we compare the total loss by friction, E h + E v , and E s . We assume that the energy that would be loaded elastically in absence of modulation, E s , is entirely dissipated by friction, which is somehow equivalent to the application of a principle of maximum energy release rate commonly used in fracture theory [13]. We thus write E h + E v = E s , which is reasonable for a stiff system (G µP ). We thus have a first relation between θ and λ. Then, we obtain the onset from the value of λ that minimizes θ. We get that, at the onset,

θλ 2 = 2 h [2 d b + δ∆µ/µ] together with: λ h 4 + λ h 3 + λ h 2 = α µ ∆µ µ δ h 1 + 2 σ s ρgh . (4)
where we took into account that d b /δ 1. In order to contrast predictions from our theoretical arguments with the experimental results, we adjust the wavelength measured with glass beads for several cohesions using αδ∆µ/µ 2 as single fitting parameter (Figs. 2 and3) and report the results in Table I. The importance of cohesion expressed trough the ratio σ s /ρgh 0 is shown in Table I, where the typical thickness of the layer is h 0 = 5 mm. The above analysis, along with the obtained values of σ s /ρgh 0 ≈ 1, indicates that the characteristics of the structures observed with small particles are strongly influenced by the cohesion (Fig. 2). On the contrary, for larger particles, σ s /(ρgh 0 ) 1, and null influence of the cohesion is expected in this limit. According to Eq. ( 4), all data presented in Fig. 3 collapse together within the errors bars. Finally, we report in Table II a summary of all the experimental parameters. The typical size l R of the surface asperities for glass and brass particles are obtained by atomic force measurements as discussed in [11]. The values of the E/J for glass and brass are from the literature. We assume a value of frictional coefficient µ = 0.5 and ∆µ = 0.1 in all cases and estimate the characteristic length δ. We observe that δ is of the order of a few micrometers, thus a fraction of the grain size. Even if no clear dependence of δ on the experimental parameters can be deduced from our measurements, the order of magnitude is in agreement with the physical meaning of this parameter, i.e. the relative displacement of two granular surfaces with respect to one another for the system to dilate, and thus to reach the constant value µ of the friction coefficient. Moreover, we note that such distance δ is compatible with a typical shear deformation θ at the onset of the order of δ/d, thus of a few percents as observed experimentally.

V. CONCLUSIONS

In conclusion, we have evidenced that sheardeformation of a granular material strongly depends on the cohesion. For σ s /P 1, the material breaks in a series of parallel cracks. In contrast, for σ s /P 1, the granular layer deforms without fracturing, but a square pattern consisting in a smooth, periodic, modulation of the layer thickness is observed. At low cohesion the weakening of the liquid bridges, which is at play at large cohesion, cannot be invoked anymore. We show that the decrease of the friction, associated with the shear-induced dilation of the material, can explain the instability observed in dry materials. Our model recovers the correct dependence of the characteristic length of the pattern on the layer thickness.

FIG. 1 .

 1 FIG. 1. (Color online)-Sketch of the experimental setup -The 4 arms of a latex membrane are driven by computercontrolled actuators. At center, arrows show the displacements of the membrane. Pure shear strain is achieved to better than 1% over a surface area of 50 cm 2 . Inset: typical modulation of shear strain at the free surface [Brass beads, d = (75 -106) µm, h = 3 mm, θ = 0.14 and RH = 39 %].

FIG. 2 .

 2 FIG. 2. (Color online) a) Pattern structure for various h and humidity. RH (σs): 39 % (17 Pa), 74 % (37 Pa), 89 % (144 Pa), scale bar = 1 cm. b) Pattern wavelength, λ, for increasing h at various cohesion (relative humidity, RH ). The lines are from Eq. (4) with parameters reported in Table I [Glass beads, d = (0 -45) µm, ∆θ = 1.5 × 10 -3 ].

FIG. 3 .

 3 FIG. 3. (Color online) a) Pattern structure for various h and humidity RH (σs): 39 % (0.5Pa), 74 % (0.8Pa), 89 %(1.9Pa), scale bar = 1 cm. b) Wavelength λ for increasing h at various cohesion. The lines are from Eq. (4) with parameters reported in Table I [Glass beads, diameter (150-200) µm, ∆θ = 1.5 × 10 -3 ].

3 FIG. 4 .

 34 FIG. 4. (Color online) Surface shear, xy and vorticity, ωz, fields for two distinct values of thickness, for θ = 0.083 [Glass beads, diameter (150-200) µm, ∆θ = 1.5 × 10 -3 , RH = 39 %].

FIG. 5 .

 5 FIG. 5. (Color online)-Amplitude < ∆d > of the modulation vs. shear-strain θ -a) RH = 39% and b) RH = 74% [Glass beads, diameter (0-45) µm, ∆θ = 1.5 × 10 -3 ].

FIG. 6 .

 6 FIG. 6. (Color online)-Vorticity ωz for brass particlesa) d = (75 -106) µm [h = 5 mm, θ = 0.14 and RH = 39 %]. b) Wavelength λ vs. thickness h, ( ) : d = (75 -106) µm and (•) : d = (212 -300) µm. The lines are from Eq. (4) with parameters reported in Table II [Brass beads, ∆θ = 3.6 × 10 -3 , θ = 0.14 and RH = 39 %].

TABLE I .

 I Parameter αδ∆µ/µ 2 vs. σs for glass particles.

	Glass			(0-45) µm		(150-200) µm
	σs (Pa)	17	37	144	0.5	0.8	1.9
	σs/ρgh0	0.17 0.37 1.44 0.005 0.008 0.019
	αδ∆µ µ 2	(mm) 2.85	5.0	-	8.3	14.1 13.6

TABLE II .

 II Summary of the experimental parameters.

			Glass	Brass
	d (µm)	0 -45 150 -200 75 -106 212 -300
	lR (µm)	0.1	0.07	1	1
	E/J	25 -40 25 -40 250 -300 250 -300
	α	1700	7500	3000	8500
	δ (µm)	∼ 5	∼ 4.4	∼ 8	∼ 11

The authors acknowledge the support from LIA-MSD, CNRS-France, and Conicyt-Chile through the program Fondecyt Anillo ACT-1410. F.M. acknowledges "Convenio de Desempeño" N 0 1555 UdeSantiago and the invitation as professor by the École Normale Supérieure de Lyon (ENS-Lyon) and H. A. the support of Fondecyt postdoctoral program, under grant 3160341.