Antonin Chambolle

Sergio Conti

Gilles Francfort

Korn-Poincaré inequalities for functions with a small jump set

come

Introduction

Functions of bounded deformation have been introduced to study plasticity, damage and fracture models in a geometrically linear setting [Suq78, Tem83, TS81, AG80, KT83]. The space BD(Ω) is the set of functions u ∈ L 1 (Ω; R n ) such that the symmetric part of the distributional gradient Eu = (Du+Du T )/2 is a bounded measure. They share many properties with the functions of bounded variation, for example the strain can be decomposed in a part absolutely continuous with respect to the Lebesgue measure L n , a jump part and an intermediate, so-called Cantor part,

Eu = e(u)L n + [u] ⊗ n + n ⊗ [u] 2 H n-1 J u + E c u , (1.1)
where the jump set J u is a (n -1)-rectifiable subset of Ω, [u] : J u → R n denotes the jump of u, and n the normal to J u , see [START_REF] Ambrosio | Fine properties of functions with bounded deformation[END_REF]. The three terms in (1.1) are mutually orthogonal and have a clearly distinct physical interpretation: e(u) represents the regular strain, [u] the crack opening or the plastic slips, E c u some notion of diffuse damage. The modeling of fracture in linear elasticity [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF] focuses on the interplay between the regular and the jump part, and is normally restricted to the special functions of bounded deformation SBD, defined as those u ∈ BD for which E c u = 0, see also [START_REF] Chambolle | A density result in two-dimensional linearized elasticity, and applications[END_REF][START_REF] Schmidt | Eigenfracture: an eigendeformation approach to variational fracture[END_REF][START_REF] Focardi | Asymptotic analysis of Ambrosio-Tortorelli energies in linearized elasticity[END_REF][START_REF] Iurlano | A density result for GSBD and its application to the approximation of brittle fracture energies[END_REF]. In fracture models it is natural to relate e(u) L 2 to an elastic energy and the total area of the crack H n-1 (J u ) to the fracture energy, considering functionals of the type e(u) 2 L 2 + H n-1 (J u ) which constitute the vectorial counterpart to the Mumford-Shah functional from image segmentation [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]. This leads to the study of the space SBD p (Ω), which is defined as the set of u ∈ BD(Ω) such that E c u = 0, e(u) ∈ L p and H n-1 (J u ) < ∞, see [START_REF] Bellettini | Compactness and lower semicontinuity properties in SBD(Ω)[END_REF][START_REF] Chambolle | An approximation result for special functions with bounded deformation[END_REF][START_REF] Chambolle | Addendum to: "An approximation result for special functions with bounded deformation[END_REF]. The key additional difficulty with respect to scalar models based on functions of bounded variation is the lack of control on the skew-symmetric part of the distributional gradient, Du -Du T .

Korn's inequality is one key ingredient in the study of linear elasticity and of functions of bounded deformation. In its standard version it states that, if u ∈ W 1,p (Ω; R n ) with p ∈ (1, ∞) and Ω ⊂ R n a bounded, connected Lipschitz set, then there is A ∈ R n×n such that

Du -A L p (Ω) ≤ c(p, Ω) Du + Du T L p (Ω) , (1.2)
see for example [START_REF] Rešetnjak | Estimates for certain differential operators with finite-dimensional kernel[END_REF][START_REF] Ting | Generalized Korn's inequalities[END_REF] or [START_REF] Nitsche | On Korn's second inequality[END_REF] for a proof. This has been generalized to many different settings, including the geometrically nonlinear case [START_REF] Friesecke | A theorem on geometric rigidity and the derivation of nonlinear plate theory from threedimensional elasticity[END_REF], mixed growth [START_REF] Conti | Korn's second inequality and geometric rigidity with mixed growth conditions[END_REF] and incompatible strain fields which contain dislocations [START_REF] Neff | Maxwell meets Korn: a new coercive inequality for tensor fields in R N ×N with squareintegrable exterior derivative[END_REF][START_REF] Müller | Geometric rigidity for incompatible fields and an application to strain-gradient plasticity[END_REF]. Korn's inequality (1.2) does not hold for p = 1, and indeed one can construct functions such that Du + Du T ∈ L 1 but Du ∈ L 1 , see [START_REF] Ornstein | A non-inequality for differential operators in the L 1 -norm[END_REF][START_REF] Conti | A new approach to counterexamples to L 1 estimates: Korn's inequality, geometric rigidity, and regularity for gradients of separately convex functions[END_REF][START_REF] Kirchheim | Automatic convexity of rank-1 convex functions[END_REF]. Correspondingly, BD is not a subset of BV , in the sense that for any open set Ω ⊂ R n there are functions u ∈ L 1 (Ω; R n ) such that Du + Du T is a bounded measure, but Du is not.

A combination of Poincaré's and Sobolev inequalities with Korn's inequality (1.2) leads to estimates on u(x) -Axb in L q , q ≥ 1. Although Korn's inequality fails for p = 1, the combination with the inequalities in which the derivative is lost is still true. Indeed, for u ∈ BD(Ω) there are A ∈ R n×n and b ∈ R n such that u(x) -(Ax + b) L 1 * (Ω) ≤ c(Ω)|Eu|(Ω) ,

(1.3) see [START_REF] Kohn | New integral estimates for deformations in terms of their nonlinear strains[END_REF][START_REF] Temam | Problèmes mathématiques en plasticité[END_REF]. Here 1 * = n/(n -1) is the Sobolev conjugate exponent to 1, and Ω ⊂ R n is again a bounded, connected Lipschitz set. For functions in SBD p one may expect a stronger result to hold, since a stronger control of the gradient is present. From a physical viewpoint, we consider a sample which contains small fractures, as illustrated in Figure 1. If the total (n -1)-dimensional area of the fracture is small, most of the material still belongs to one large piece, whose deformation is controlled only by the regular strain e(u). Small subsets can however be completely detached, therefore we cannot expect an estimate in the entire set Ω, but only in a subset Ω \ ω. The "holes" ω are the parts of Ω which are effectively separated from the rest by J u . For their volume we obtain first the simpler estimate L n (ω) ≤ cH n-1 (J u )diam (Ω) in Section 2 and then the optimal isoperimetric estimate L n (ω) ≤ c(H n-1 (J u )) n/(n-1) in Section 4. Further, a corresponding trace estimate can be obtained, see Section 4. Our bound on u in L q with q = pn/(n -1) only reaches the optimal exponent of the full Korn-Sobolev inequality (1.3) if p = 1.

The key idea of the proof is to use the fundamental theorem of calculus along lines which do not intersect the jump set to estimate the variation of u. This would be immediate in the BV setting, in which one fully controls the gradient. In a BD setting one only obtains control of the longitudinal component of u. Therefore one needs to consider lines with many different orientations, making sure to choose them so that they do not intersect the jump set, a strategy that was used for proving density in [START_REF] Chambolle | An approximation result for special functions with bounded deformation[END_REF][START_REF] Chambolle | Addendum to: "An approximation result for special functions with bounded deformation[END_REF] and for proving rigidity in one and multiwell settings for example in [START_REF] Kohn | New integral estimates for deformations in terms of their nonlinear strains[END_REF][START_REF] Dolzmann | Microstructures with finite surface energy: the two-well problem[END_REF][START_REF] Conti | Rigidity and Gamma convergence for solid-solid phase transitions with SO(2) invariance[END_REF]. More details are explained in the introduction to Section 2.

Our main result is the following.

Theorem 1. Let Ω ⊂ R n be a connected bounded Lipschitz set, p ∈ [1, ∞).
There is a constant c which depends only on p and Ω such that for any u ∈ SBD(Ω), one can find Borel sets ω ⊂ Ω and ω Γ ⊂ ∂Ω with

L n (ω) ≤ c(H n-1 (J u )) n/(n-1) and H n-1 (ω Γ ) ≤ cH n-1 (J u ) (1.4)
and an affine function a : R n → R n with Da + Da T = 0 such that, with q = pn/(n -1),

Ω\ω |u -a| q dx ≤ c Ω |e(u)| p dx n/(n-1)
(1.5)

and ∂Ω\ω Γ |T u -a| p dH n-1 ≤ c Ω |e(u)| p dx . (1.6)
Here e(u) denotes the part of the strain Eu = (Du + Du T )/2 which is absolutely continuous with respect to L n and T u denotes the trace of u on ∂Ω.

From the statement it is clear that the assertion is only relevant if u ∈ SBD p for some p ≥ 1; the case p = 1 is included.

The original motivation for this work was in the study of a fracture model with a non-interpenetration constraint, which will be discussed elsewhere [CCF]. A different application, where the trace estimate is specifically used, arises in the study of pattern formation in delaminated thin films, see [BCM].

A natural question is whether a similar estimate for ∇u can be derived. We do not know if this is true and, in particular, we do now know if SBD p , p > 1, is a subset of BV .

The rest of this paper contains the proof of Theorem 1. In Section 2 we construct the affine map a with suboptimal estimates on both u and |ω|. In Section 3 we build upon this result to prove higher integrability of u in the interior. In Section 4 we refine the argument of Section 3 to prove higher integrability up to the boundary and the trace estimate; the proof of Theorem 1 then follows easily by covering.

Notation. As stated above, for Ω ⊂ R n open BD(Ω) is the set of functions u ∈ L 1 such that Eu = (Du+Du T )/2 is a bounded measure. Then Eu obeys (1.1), and SBD is the subset of BD such that E c u = 0. We denote by e(u) the density of Eu with respect to L n , as in (1.1).

We shall use the results on slicing of BD functions from [START_REF] Ambrosio | Fine properties of functions with bounded deformation[END_REF]. In particular, let u ∈ SBD, ξ ∈ S n-1 . We set u ξ y (t) = ξ • u(y + tξ). By the structure theorem [ACD97, Th. 4.5], for almost every y one has: u ξ y ∈ SBV , du ξ y /dt = ξe(u)ξ L 1 -a.e., and the section of the jump set is the jump set of u ξ y . In particular we have that, for a.e.

x such that x + [0, 1]ξ ⊂ Ω \ J u , ξ • u(x) -ξ • u(x + ξ) = 1 0 ξ • e(u)(x + tξ)ξdt . (1.7)
Functions in BD have an L 1 trace on the boundary, see [Suq79, TS81, Bab13], which we denote by T u.

Throughout the paper c denotes a generic constant, which may change from line to line, and which, unless otherwise stated, depends only on the spatial dimension n and on the integrability exponent p.

Acknowledgements The authors gratefully acknowledge the hospitality of the Mathematisches Forschungsinstitüt Oberwolfach (MFO) where this research was initiated. A. Chambolle is partially funded by the Agence Nationale de la Recherche, project ANR-12-BS01-0014-01 "GEOMETRYA".

The local estimate

We start proving the Korn-Poincaré estimate in a square, with an exceptional set whose volume is proportional to H n-1 (J u ). This is the part of the proof where we construct the affine map a. Before stating our result we give a sketch of the main ideas.

The key strategy is to use the fundamental theorem of calculus along segments which do not intersect the jump set, in order to relate the values of u at different points. Since we only control the strain and not the full gradient we need to focus on the longitudinal component, and with enough regularity we compute

d dt ξ • u(x + tξ) = ξ • Du(x + tξ)ξ = ξ • Du + Du T 2 (x + tξ)ξ . (2.1)
We shall use this to relate the values of u in n + 1 points, called {z 0 , . . . , z n }, which form the vertices of an n-dimensional simplex, see Figure 2. Precisely, integrating (2.1) gives a control on the longitudinal variations,

|(z i -z j ) • (u(z i ) -u(z j )| ≤ |z i -z j | [z i ,z j ] |e(u)|dH 1
whenever the restriction of u to the line through z i and z j belongs to W 1,1 . We shall use the n + 1 values u(z i ) to define a linear map ã : R n → R n with ã(z i ) = u(z i ). For this map as well, the longitudinal variations are related to the components of the symmetrized gradient. The longitudinal variations are those of u, hence they remain small. Therefore, since ã is affine, the symmetric part of the gradient of ã is also bounded by e(u) L 1 . We can then replace ã by its projection on affine functions with skew-symmetric gradients.

z 0 z 1 z 2 y
Figure 2: Sketch of the construction in Proposition 2. We first obtain rigidity on the vertices z i = z 0 +t * e i of a simplex, and then show that the same rigidity is inherited by the other points y ∈ Q.

This strategy works under three assumptions: that the segments [z i , z j ] do not intersect the jump set, that we can use the fundamental theorem of calculus along these segments, and that the L 1 norm of e(u) on the segments can be estimated by the L 1 norm on the entire set. The first and the third are true "on average", if we average over sets of parallel segments first, and therefore for many possible choices of the z i . The second one is true by the structure theorem for almost every choice.

It remains to show that u is close to a on a large part of the square. To do this we consider a generic point y and repeat the above argument on the segments [y, z i ]. This leads to an estimate on the n + 1 longitudinal components of (u-a)(y) in terms of the L 1 norm of e(u) along those segments. A simple argument from convex analysis shows that these components are linearly independent in a uniform way, permitting to estimate |u -a|(y) from them. Finally we need to make sure that we can use the fundamental theorem of calculus for many points y and that the integral of e(u) along the segments [y, z i ] can be estimated by the L 1 norm of e(u). Both require additional conditions on the choices of the point z 0 , see Figure 3. The final choice will have to fulfill the four conditions given in (2.8), (2.9), (2.12), and (2.17). (2.

Proposition 2. Let u ∈ SBD(Q r ), Q r = (-r, r) n , p ∈ [1, ∞). Then there are ω ⊂ Q r with L n (ω) ≤ crH n-1 (J u ) (2.2) J u (a) J u y (b)
3)

The constant depends only on n and p.

Proof. By rescaling we can assume r = 1; we write for brevity

Q = Q 1 .
Instead of finding a with e(a) = 0 it suffices to construct an affine function a with |e(a)| ≤ c e(u) L 1 , which will be proven in (2.11), and then to replace a by its projection on infinitesimal rotations. Further, we can assume

H n-1 (J u ) ≤ 1 32n 3 .
(2.4) Indeed, if not, one can simply take ω = Q and a = 0 with c = 32 • n 3 . We subdivide the proof in several steps.

Step 1: Fundamental theorem of calculus along segments.

We consider for (x, ξ, t)

∈ Q × S n-1 × R the relation ξ • (u(x + tξ) -u(x)) = t 1 0 ξ • e(u)(x + stξ)ξ ds . (2.5) If u ∈ C 1 (Q; R n ) this would hold for all (x, ξ, t) with x+tξ ∈ Q.
Here we need to pay attention to both to the jump set and to the possible lack of pointwise values and differentiability on the appropriate null sets. To do this we define the characteristic function of the exceptional set T :

R n × S n-1 × R → R by T (x, ξ, t) = 1 if x ∈ Q , x + tξ ∈ Q and (2.5) does not hold, 0 otherwise. (2.6)
Both sides in (2.5) are measurable, hence T is measurable. We inserted the condition x, x + tξ ∈ Q since we are only interested in points inside Q. The condition "(2.5) does not hold" includes the case in which the integral is not defined.

By the structure theorem [ACD97, Th. 4.5], for any fixed ξ the following holds: for almost every x ∈ Q the section t → ξ • u(x + tξ) is in SBV , its jump set coincides almost everywhere with the section of the jump set of u, and its derivative is ξ • e(u)ξ. This means that for L n -almost every x ∈ Q one of these two options holds: either x + Rξ intersects J u or (2.5) holds for L 1 -almost every t in the relevant interval. Since, for any ξ ∈ S n-1 , the set

ω ξ = {x ∈ Q : x + Rξ ∩ J u = ∅} has measure bounded by |ω ξ | ≤ H n-1 (J u )diam (Q) = 2 √ nH n-1 (J u ), we obtain the estimate R n R T (x, ξ, t)dtdx ≤ |ω ξ |diam (Q) ≤ 4nH n-1 (J u )
for any ξ ∈ S n-1 .

(2.7) We need now to choose the point z 0 as well as n other reference points, obtained from the point z 0 by shifts along the canonical basis, and we need to control also their interactions, see Figure 2. Therefore we define, for

z 0 ∈ q = (-1, 0) n ⊂ Q and t ∈ (1/2, 1), G(z 0 , t) = n i=1 T (z 0 , e i , t) + n i=1 n j=i+1 T (z 0 + te i , e j -e i √ 2 , t √ 2) .
Integrating and using (2.7) on each term leads to

1 1/2 q G(z 0 , t)dz 0 dt ≤ R R n Gdz 0 dt ≤4n 2 H n-1 (J u ) + n(n -1) 2 4n √ 2 H n-1 (J u ) ≤ 4n 3 H n-1 (J u ) .
We choose t * ∈ (1/2, 1) at a set value for the rest of the proof, such that q G(z 0 , t * )dz 0 ≤ 8n 3 H n-1 (J u ). Therefore, (2.4) implies that we can have G(z 0 , t * ) = 0 only for at most one-quarter of the values of z 0 ∈ q, and therefore the statement (2.5) holds on all edges [z 0 , z 0 + t * e i ], i = 1, . . . , n and [z 0 + t * e i , z 0 + t * e j ], 1 ≤ i < j ≤ n (2.8) is true for three-quarters of the values of z 0 ∈ q. The choice of z 0 among these possibilities will occur later.

Step 2: Estimate for e(u) along the edges of the simplex, construction of a.

For t * as above and any z 0 ∈ q we set z i = z 0 + t * e i ∈ Q, i = 1, . . . , n, and define

F (z 0 ) = n i,j=0 [z i ,z j ] |e(u)|dH 1 .
We integrate in z 0 , change variables and the order of integration to obtain (setting e 0 = 0)

q F (z 0 )dz 0 = n i,j=0 q t * 0 |e(u)|(z 0 + t * e i + s(e j -e i ))|e i -e j |dsdz 0 ≤ √ 2 n i,j=0 t * 0 Q |e(u)|(x)dxds ≤ √ 2(n + 1) 2 e(u) L 1 (Q) .
Therefore for three-quarters of the choices of z 0 ∈ q one has

F (z 0 ) ≤ 4 √ 2(n + 1) 2 e(u) L 1 (Q) .
(2.9) This is the second requirement on z 0 , and we still have at least half the volume of q at our disposal to enforce additional conditions later.

We consider any choice of z 0 such that both properties (2.8) and (2.9) hold. Then (2.5) gives

|(e i -e j ) • (u(z i ) -u(z j ))| ≤ c e(u) L 1 (Q) , i, j = 0, . . . , n .
(2.10)

We define the affine map a : R n → R n by setting a(z i ) = u(z i ), i = 0, . . . , n. From (2.10) we obtain

t * |(Da) ii | = |e i • (a(z 0 + t * e i ) -a(z 0 ))| ≤ c e(u) L 1 (Q)
as well as

t * |(Da) ii + (Da) jj -(Da) ij -(Da) ji | = |(e i -e j ) • (a(z 0 + t * e i ) -a(z 0 + t * e j ))| ≤ c e(u) L 1 (Q) ,
which together yield

|Da + Da T | ≤ c e(u) L 1 (Q 1 ) .
(2.11)

We define w(x) = u(x)a(x), so that w(z i ) = 0 for i = 0, . . . , n.

Step 3: Definition of the exceptional set ω.

In order to prove a bound on |w|(y) for many points y ∈ Q we need to ensure validity of the relation (2.5) for the segments [z i , y]. We first observe that for three-quarters of the choices of z 0 ∈ q one has

S n-1 R T (z i , ξ, t) dt dH n-1 (ξ) ≤ cH n-1 (J u ) (2.12)
for all i, with a constant c which depends only on n. To see this, it suffices to integrate (2.7) over ξ ∈ S n-1 and swap the order of integation. This is the third condition on z 0 , and we still have at least one-quarter of q at our disposal, which we shall use later to enforce a fourth condition.

For every i = 0, . . . , n we define the set

ω (i) = {y ∈ Q : y = z i + tξ with T (z i , ξ, t) = 1} ,
which by (2.12) it has measure bounded by

|ω (i) | = S n-1 ∞ 0 t n-1 T (z i , ξ, t) dt dH n-1 (ξ) ≤ cH n-1 (J u ) .
We define

ω = n i=0 ω (i) (2.13)
and observe that (2.5) holds for any segment [z i , y] with y ∈ Q \ ω and i ∈ {0, . . . , n}. Since w differs from u only by the affine function a, (2.5) holds also for w.

Step 4: Estimate for w on Q \ ω.

It remains to ensure integrability of the strain along the relevant segments. We write, using (2.5) on each segment [y, z i ],

n i=0 |(y -z i ) • (w(y) -w(z i ))| ≤ n i=0 |y -z i | [y,z i ] |e(w)|dH 1 . (2.14)
Since w(z i ) = 0, if the vectors yz i were the basis vectors we would immediately obtain an estimate on |w|(y). To see that the same is true for the given estimate, we need to show that these n + 1 vectors are "uniformly linearly independent", in a sense we now make clear. We first observe that the convex hull C of the n + 2 points {y, z 0 , . . . , z n } has volume larger than t n * /n! (since this is the volume of the convex hull of the {z i } alone). We now claim that C coincides with the union of the convex hull of the n + 1 sets {y, z 0 , . . . , z n } \ {z i }. To see this, pick any p ∈ C. By Caratheodory's theorem, p is a convex combination of n + 1 of the points. If y is one of them, we are done. Otherwise, p ∈ K = conv {z 0 , . . . , z n }. Let q be the last point of the ray y + [0, ∞)(py) contained in K. Since q is on the boundary of K, it is a convex combination of n of the z i , therefore p is a convex combination of y and n of the z i , concluding the proof of the claim.

Therefore for any y there are n of the points {z 0 , . . . , z n } such that the simplex they generate has volume at least t n * /(n + 1)!. This means, that the n corresponding vectors yz i have the property that the determinant of the n×n matrix of which they are the columns is at least t n * /(n+1) ≥ 2 -n /(n+1). This n × n matrix A has coefficients a i,j ∈ [-2, 2] for all i, j, so that its comatrix is bounded by some constant C(n) depending only on n. It follows

|A -1 | ≤ C(n)/| det A| and therefore | det A| C(n) |x| ≤ |Ax| for all x ∈ R n .
Consequently, bounding the left-hand side of (2.14) from below and using that w(z i ) = 0, we obtain

|w(y)| ≤ c n i=0 [y,z i ] |e(w)|dH 1 for all y ∈ Q \ ω ,
and with Hölder's inequality

|w| p (y) ≤ c n i=0 [y,z i ] |e(w)| p dH 1 for all y ∈ Q \ ω .
(2.15)

It remains to integrate in y. To estimate the contribution of the i-th term in (2.15) we set g = |e(w)| p χ Q and change variables from

y to ỹ = y -z 0 ∈ Q 2 = (-2, 2) n , Q\ω [y,z i ] |e(w)| p dH 1 dy ≤ Q 2 [z 0 +ỹ,z 0 +e i ] g(x)dH 1 (x)dỹ = Q 2 [ỹ,e i ]
g(z 0 + t)dH 1 (t)dỹ .

We sum over i and set

H(z 0 ) = n i=0 Q 2 [ỹ,e i ]
g(z 0 + t)dH 1 (t)dỹ , so that Q\ω |w| p (y)dy ≤ cH(z 0 ) .

(2.16)

Swapping the order of integration, we estimate as usual

q H(z 0 )dz 0 ≤ n i=0 Q 2 [e i ,ỹ] g L 1 (R n ) dH 1 (t)dỹ ≤ c e(w) p L p (Q 1 )
where c depends only on n. Therefore there are many choices of z 0 (more than three-quarters) such that

H(z 0 ) ≤ 5c e(w) p L p (Q 1 ) .
(2.17) This is the fourth and last condition we impose on z 0 . At this point we can fix z 0 such that the four conditions (2.8), (2.9), (2.12), and (2.17) are fulfilled.

Recalling that (2.11) gives e(w) L p (Q 1 \ω) ≤ c e(u) L p (Q 1 ) , (2.16) and (2.17 and conclude the proof.

Higher integrability

In this Section we show that control of e(u) in L p gives an estimate for ua in a space better than L p . However, we only get the estimate in the optimal space L p * if p = 1. The basic strategy is to use rigidity in finitely many directions ξ, and then to combine estimates along different directions in a way which is similar to the one used in proving the Sobolev embedding of W 1,1 into L (n-1)/n , see Lemma 4 below. The standard argument for u ∈ W 1,p , based upon considering f = |u| p ∈ W 1,1 which obeys |Df | ≤ p|u| p-1 |Du|, does not apply here since we only control u on part of the domain. (i) There exist a set ω ⊂ Q ′ and an affine function a : R n → R n with e(a) = 0 such that

Proposition 3. Let Q = (-r, r) n , Q ′ = (-r/2, r/2) n , u ∈ SBD(Q), p ∈ [1, ∞). R x ξ x R y ξ y Q Q ′
L n (ω) ≤ crH n-1 (J u ) (3.1)
and

Q ′ \ω |u -a| np/(n-1) dx ≤ cr n(p-1)/(n-1) Q |e(u)| p dx n/(n-1)
.

(ii) If additionally p > 1 then there is q > 0 (depending on p and n) such that, for a given mollifier ϕ r ∈ C ∞ c (B r/4 ), ϕ r (x) = r -n ϕ 1 (x/r), the function v = uχ Q ′ \ω + aχ ω obeys

Q ′′ |e(v * ϕ r ) -e(u) * ϕ r | p dx ≤ c H n-1 (J u ) r n-1 q Q |e(u)| p dx ,
where Q ′′ = (-r/4, r/4) n .

The constant in (i) depends only on p and n, the one in (ii) also on ϕ 1 .

Proof. By Proposition 2 applied to Q = (-r, r) n there exist a set ω 0 ⊂ Q with |ω 0 | ≤ crH n-1 (J u ) and an affine function a with e(a) = 0 such that

Q\ω 0 |u -a| p dx ≤ cr p Q |e(u)| p dx . (3.3)
We can assume without loss of generality that a = 0 (otherwise we replace u by ua) and r = 1 (by scaling). For ξ ∈ S n-1 and x ∈ Q ′ we define the ray

R x ξ = (x + Rξ) ∩ Q ,
see Figure 4. Then, R x ξ is a segment and H 1 (R x ξ ) ≥ 1. For any direction ξ ∈ S n-1 , the structure theorem [ACD97, Th. 4.5] states that for almost every z for which R z ξ ∩ J u = ∅ one has, in a way which closely corresponds to (2.5) in the proof of Proposition 2,

ξ • (u(y) -u(x)) = 1 0 ξ • e(u)(x + s(y -x))(y -x)ds for H 1 -a.e. x, y ∈ R z ξ . (3.4)
For ξ ∈ S n-1 we define the "shadow" of the jump set

ω ξ = {z ∈ Q ′ : (3.4) does not hold} .
This corresponds, up to a null set, to the points z such that a line through z parallel to ξ intersects the jump set, hence to the points in which we cannot use the estimate in direction ξ, and it obeys H n-1 (Π ξ ω ξ ) ≤ H n-1 (J u ), where Π ξ : R n → R n denotes the orthogonal projection onto the (n -1)dimensional space ξ ⊥ . We remark that this set is cylindrical, in the sense that

ω ξ = (ω ξ + ξR) ∩ Q ′ .
Given a direction ξ ∈ S n-1 we define the exceptional set

ω * ξ = ω ξ ∪ {z ∈ Q ′ : H 1 (R z ξ ∩ ω 0 ) ≥ 1/2},
where ω 0 is the set entering (3.3), which was obtained from Proposition 2. It is easy to see that

H n-1 (Π ξ ω * ξ ) ≤ cH n-1 (J u ). For almost all x ∈ Q ′ \ ω * ξ , (3.4) gives |ξ • u(x)| ≤ |ξ • u(y)| + R x ξ |e(u)|dH 1 for a.e. y ∈ R x ξ .
Averaging over all y ∈ R x ξ \ ω 0 gives

|ξ • u|(x) ≤ 1 H 1 (R x ξ \ ω 0 ) R x ξ \ω 0 |u|dH 1 + R x ξ |e(u)|dH 1 and |ξ • u|(x) ≤ 2 x+Rξ f dH 1 for a.e. x ∈ Q ′ \ ω * ξ , (3.5) where f = |e(u)|χ Q + |u|χ Q\ω 0 .
We observe that by (3.3) the function f obeys

f L p (R n ) ≤ c e(u) L p (Q) . (3.6)
We choose a set S containing n 2 vectors of S n-1 with the property that any n of them are linearly independent. To show that they exist it suffices to choose them iteratively, at any step only a null set needs to be avoided. These vectors are universal and will influence the constants in the estimate. At variance with Proposition 2 we shall use here only finitely many directions.

For any set V ⊂ S of n vectors (which are automatically linearly independent) one has

|u|(x) ≤ c V ξ∈V |u • ξ|(x) ,
with a constant which depends only on V . We define

ω = ξ∈S ω * ξ which obeys ω ⊂ Q ′ , |ω| ≤ cH n-1 (J u
) and therefore (3.1). Recalling (3.5),

|u|(x) ≤ c V ξ∈V x+Rξ f dH 1 for a.e. x ∈ Q ′ \ ω
and by Hölder's inequality

|u| p (x) ≤ c V ξ∈V x+Rξ f p dH 1 for a.e. x ∈ Q ′ \ ω. (3.7)
At this point we use the estimate (3.7) with several different sets of vectors. Let V 1 , . . . , V n be n disjoint sets of n vectors in S. Then (3.7) holds for each of them. We multiply these n equations and obtain for almost every

x ∈ Q ′ \ ω |u| np (x) ≤ n j=1 c V j ξ∈V j x+Rξ f p dH 1 for a.e. x ∈ Q ′ \ ω. (3.8)
Since there are finitely many choices of V ⊂ S the constant is universal. We swap the sum with the product and obtain

|u| np (x) ≤ c S W ⊂S,#W =n ξ∈W x+Rξ f p dH 1 for a.e. x ∈ Q ′ \ ω,
where by the definition of S the vectors in any of the sets W = {ξ 1 , . . . , ξ n } ⊂ S are linearly independent (the relevant sets W are those which contain one vector from each of the V j ). By Lemma 4 (see below) applied to each term in the sum and f p we obtain

uχ Q ′ \ω L np/(n-1) ≤ c f L p .
Recalling (3.6) the proof of (3.2) and therefore of part (i) is concluded. We now turn to assertion (ii) concerning the function v = uχ Q ′ \ω . We define, with

ϕ = ϕ 1 ∈ C ∞ c (B 1/4 ) as in the statement, a function h ∈ C ∞ (Q ′′ , R n×n sym ) by h = e(v * ϕ) -e(u) * ϕ . (3.9)
To conclude the proof it suffices to bound ξ • hξ in L p (Q ′′ ) for a finite set of ξ. For example, all ξ ∈ T = {e i } i=1,...,n ∪ {(e i + e j )/ √ 2} i,j=1...n will do, since A → ξ∈T |ξ • Aξ| is a norm on R n×n sym . We fix ξ ∈ S n-1 and define, for η ∈ (0, 1/2) chosen below,

ω * * ξ = ω ξ ∪ {z ∈ Q ′ : H 1 (R z ξ ∩ (ω ∪ ω 0 )) ≥ η}, which obeys ω * ξ ⊂ ω * * ξ and |ω * * ξ | ≤ c H n-1 (J u ) η . (3.10) We pick z ∈ Q ′ \ ω * * ξ . Then z ∈ ω * ξ and, since the latter set is cylindrical, R z ξ ∩ ω * ξ = ∅. Therefore (3.5) gives |ξ • u|(x) ≤ 2 R z ξ f dH 1 for a.e. x ∈ Q ′ ∩ R z ξ ,
and, with Hölder's inequality,

|ξ • u| p (x) ≤ c R z ξ f p dH 1 for a.e. x ∈ Q ′ ∩ R z ξ . Integration over all x ∈ ω ∩ R z ξ ⊂ Q ′ ∩ R z ξ gives, since H 1 (ω ∩ R z ξ ) < η, ω∩R z ξ |ξ • u| p dH 1 ≤ cη R z ξ f p dH 1 .
Integrating then over z ∈ Q ′ \ ω * * ξ and using (3.6) yields

ω\ω * * ξ |ξ • u| p dx ≤ cη Q f p dx ≤ cη e(u) p L p (Q) . (3.11)
For the same ξ we define the longitudinal component w (ξ) : R n → R by

w (ξ) = ξ • uχ Q ′ \ω * * ξ .
We now show that for any z

∈ Q ′ the function t → w (ξ) (z + tξ) is in W 1,1 (I), with I = {t ∈ R : z + tξ ∈ Q ′ }, and
d dt w (ξ) (z + tξ) = ξ • e(u)ξχ Q ′ \ω * * ξ (z + tξ) . (3.12)
To see this, we observe that ω * * ξ is cylindrical, in the sense that

ω * * ξ = (ω * * ξ + ξR) ∩ Q ′ . If z ∈ ω * *
ξ then w (ξ) (z + tξ) = 0 for all t and (3.12) holds. If instead z ∈ ω * * ξ then χ Q ′ \ω * * ξ (z + tξ) = 1 for all t ∈ I, so that w (ξ) = ξ • u. Since z ∈ ω ξ (3.4) holds, and this implies t → ξ • u(z + tξ) ∈ W 1,1 (I) with derivative ξ • e(u)ξ(z + tξ). This proves (3.12).

In particular,

D ξ w (ξ) = ξ • e(u)ξχ Q ′ \ω * * ξ
and therefore in Q ′′ we have, recalling (3.9),

ξ • hξ = D ξ (ξ • v * ϕ) -ξ • e(u)ξ * ϕ = D ξ ((ξ • v -w (ξ) ) * ϕ) + D ξ (w (ξ) * ϕ) -ξ • e(u)ξ * ϕ = D ξ ((ξ • v -w (ξ) ) * ϕ) -ξ • (e(u)χ ω * * ξ )ξ * ϕ . (3.13)
We estimate the two terms in (3.13) separately in L p (Q ′′ ).

In the first term we use |Dϕ| ≤ c to obtain

Q ′′ D ξ (ξ • v -w (ξ) ) * ϕ p dx ≤ c Q ′ |ξ • v -w (ξ) | p dx . Since w (ξ) = ξ • v = ξ • u on Q ′ \ (ω ∪ ω * * ξ ), w (ξ) = 0 on ω * * ξ , and ξ • v = 0 on ω, Q ′ |ξ • v -w (ξ) | p dx = ω\ω * * ξ |w (ξ) | p dx + ω * * ξ \ω |ξ • v| p dx .
The last term, using Hölder's inequality, (3.10) and (3.2), is controlled by

ω * * ξ \ω |ξ • v| p ≤|ω * * ξ | 1/n Q ′ \ω |v| np/(n-1) (n-1)/n ≤c H n-1 (J u ) η 1/n e(u) p L p (Q) .
Instead, the first one is controlled by (3.11), which gives

ω\ω * * ξ |w (ξ) | p dx ≤ cη e(u) p L p (Q) .
In the second term in (3.13) we estimate

e(u)χ ω * * ξ * ϕ L p (Q ′′ ) ≤ e(u)χ ω * * ξ L 1 (Q ′ ) ϕ L p (R n ) ≤ χ ω * * ξ L p ′ (Q ′ ) e(u) L p (Q ′ ) ϕ L p (R n ) ≤ c H n-1 (J u ) η 1-1/p e(u) L p (Q ′ ) ,
where we used Hölder's inequality with p ′ defined by 1/p + 1/p ′ = 1 and (3.10).

Combining these estimates with (3.13) leads to

ξ•hξ L p (Q ′′ ) ≤ c H n-1 (J u ) η 1/(np) + η 1/p + H n-1 (J u ) η 1-1/p e(u) L p (Q) .
If 1/(np) ≤ 1 -1/p we choose η = (H n-1 (J u )) 1/(n+1) and obtain

h L p (Q ′′ ) ≤ c(H n-1 (J u )) 1/(p(n+1)) e(u) L p (Q) .
If instead 1/(np) > 1 -1/p we choose η = (H n-1 (J u )) (p-1)/p and obtain

h L p (Q ′′ ) ≤ c(H n-1 (J u )) (p-1)/p 2 e(u) L p (Q) .
Lemma 4. Let V ⊂ S n-1 be a set of n linearly independent vectors, u, f ∈ L 1 (R n ) obey

|u| n (x) ≤ n i=1 x+Rξ i f dH 1 (3.14)
for almost every x ∈ R n . Then u ∈ L n/(n-1) (R n ) and

u L n/(n-1) (R n ) ≤ c f L 1 (R n ) . (3.15)
Proof. This is part of the standard proof of Sobolev inequality, as discussed for example in [Eva98, Eq. (11-13), pages 263-264]. We remark that this part of the proof only deals with Fubini and iterated integration and only requires f ∈ L 1 (R n ).

Estimate up to the boundary

In this section we refine the argument from the previous one in several directions. We obtain a better estimate of the volume of the exceptional set, with the isoperimetric scaling, we obtain an estimate up to the boundary, and an estimate for the trace. A covering argument will then lead to the proof of Theorem 1. One important ingredient in proving a better estimate for the volume of the exceptional set is the following isoperimetric estimate, which bounds the volume of a set from the area of its codimension-1 projections.

Lemma 5. Let n ≥ 2, Ω ⊂ R n be measurable and bounded, let ξ 1 . . . , ξ n ∈ S n-1 be linearly independent vectors. Then

L n (Ω) ≤ c n i=1 (H n-1 (Π ξ i Ω)) 1/(n-1) .
(4.1)

The constant c depends only on the vectors ξ. Here Π ξ : R n → R n denotes the projection onto the n -

1-dimensional space ξ ⊥ , Π ξ = Id -ξ ⊗ ξ/|ξ| 2 .
Proof. After a linear transformation we can assume ξ i = e i . We shall prove that

L n (Ω) ≤ n i=1 H n-1 (Π e i Ω) 1/(n-1) . (4.2)
We write

ω i = Π e i Ω. If n = 2 one has χ Ω (x 1 , x 2 ) ≤ χ ω 2 (x 1 , 0)χ ω 1 (0, x 2 )
and integrating with Fubini's theorem gives

L 2 (Ω) = R 2 χ Ω dx ≤ H 1 (ω 1 )H 1 (ω 2 ) .
To prove (4.2) in higher dimension we proceed by induction on n. We define for t ∈ R

Ω t = Ω ∩ {x 1 = t} and h(t) = H n-1 (Ω t ) .
From Π e 1 Ω t ⊂ ω 1 one obtains h(t) ≤ H n-1 (ω 1 ) for all t. Projecting in the other directions we obtain

Π e i Ω t = ω i ∩ {x 1 = t} for i = 2, . . . , n .
By the inductive assumption applied to the n -1 dimensional space {x 1 = t} we obtain

h(t) ≤ n i=2 (H n-2 (ω i ∩ {x 1 = t})) 1/(n-2) .
From h(t) ≤ min{a, b} it follows that h(t) ≤ a 1/(n-1) b (n-2)/(n-1) , and therefore

h(t) ≤ (H n-1 (ω 1 )) 1/(n-1) n i=2 (H n-2 (ω i ∩ {x 1 = t})) 1/(n-1) .
With Hölder's inequality we obtain

R h(t)dt ≤ (H n-1 (ω 1 )) 1/(n-1) n i=2 R H n-2 (ω i ∩ {x 1 = t})dt 1/(n-1)
and therefore, by Fubini's theorem,

L n (Ω) = R h(t)dt ≤ n i=1 H n-1 (ω i ) 1/(n-1) .
This concludes the proof of (4.2) and therefore that of the Lemma.

We now present the main local estimate. The proof is similar to that of Proposition 3, but refined in various directions. As in Proposition 3 we estimate u(x) by controlling |ξ •u|(x) for a set of n distinct, uniformly linearly independent vectors ξ. In this case the admissible ξ are all chosen in the small set S, corresponding to a uniform "interior cone condition" for the relevant part of the boundary, so that they all point "outside", permitting to obtain the estimate on u up to the boundary and, therefore, on the trace. Using n 2 vectors instead of n gives higher integrability, as in Proposition 3. The improved estimate for the size of the exceptional set ω is obtained by noticing that a point x needs to be included in ω only if there is no choice of a set of directions ξ which permits to estimate u(x). Therefore we use first a larger (but still finite) set S of vectors, containing n 2 + n -1 elements and still contained in the cone; each subset Z ⊂ S containing n 2 vectors leads to an estimate in L np/(p-1) outside a suitable (Z-dependent) exceptional set. The set ω is then the intersection of all these exceptional sets, and its volume can be estimated with Lemma 5. Proposition 6. For r > 0, ϕ ∈ Lip ((-r, r) n-1 , [r, ∞)) with inf ϕ = r, and A : R n → R n an affine isometry we define the open sets

Ω = A{(x ′ , x n ) ∈ (-r, r) n-1 × R : -r < x n < ϕ(x ′ )} and Ω int = A{(x ′ , x n ) ∈ (-r/2, r/2) n-1 × R : 0 < x n < ϕ(x ′ )} .
Let u ∈ SBD(Ω). Then there is an affine function a : R n → R n with e(a) = 0 such that:

x n ∈ R r 0 -r x ′ ∈ R n-1 r -r r 2 -r 2 R x ξ x R y ξ y Ω Ω int e n ξ ν
ω * ξ = ω ξ ∪ {x ∈ Ω int : H 1 (R x ξ ∩ ω 0 ) ≥ 1}, which obey H n-1 (Π ξ ω ξ ) ≤ cH n-1 (J u ) and H n-1 (Π ξ ω * ξ ) ≤ cH n-1 (J u ). The same argument as in Proposition, averaging over all y ∈ Q ∩ R x ξ \ ω 0 , 3 leads to |ξ • u|(x) ≤ R x ξ |e(u)|dH 1 + 1 H 1 (Q ∩ R x ξ \ ω 0 ) Q∩R x ξ \ω 0 |u|dH 1 for a.e. x ∈ Ω int \ ω * ξ . Since H 1 (Q ∩ R x ξ \ ω) ≥ 1 we obtain |ξ • u| p (x) ≤ c R x ξ f p dH 1 for a.e. x ∈ Ω int \ ω * ξ ,
where f = |e(u)|χ Ω + |u|χ Q\ω 0 . We let S be a subset of n 2 + n -1 vectors in S with the property that any n of them are linearly independent. The choice of S depends only on L and will influence the constants in the statement. For any set V ⊂ S of n vectors we have

|u| p (x) ≤ c V ξ∈V |u • ξ| p (x) ≤ c V ξ∈V R x ξ f p dH 1 . (4.7)
For a set Z ⊂ S of n 2 vectors we define

ω Z = ξ∈Z ω * ξ .
The same argument as in Proposition 3 leads to

uχ Ω int \ω Z L np/(n-1) (R n ) ≤ c Z f L p (Ω) .
The same holds for any possible choice of Z. We set

ω = Z⊂S,#Z=n 2 ω Z = Z⊂S,#Z=n 2 ξ∈Z ω * ξ .
Since there are finitely many possible choices of Z, we conclude that uχ

Ω int \ω L np/(n-1) (R n ) ≤ Z c Z f L p (Ω) ≤ c f L p (Ω)
, which concludes the proof of (4.4).

The set ω is the set of points such that for any choice of Z they lie in at least one of the ω * ξ . Since S has n 2 + n -1 elements, if a point belongs to less than n of the ω * ξ , ξ ∈ S, then there is a choice of Z such that it does not belong to ∪ ξ∈Z ω * ξ . We conclude that any point in ω belongs to at least n of the ω * ξ , ξ ∈ S, therefore

ω ⊂ W ⊂S,#W =n ξ∈W ω * ξ = W ⊂S,#W =n ω W , ω W = ξ∈W ω * ξ .
For all ξ ∈ S we have H n-1 (Π ξ ω * ξ ) ≤ cH n-1 (J u ). For any choice of W , since the n vectors are independent Lemma 5 implies |ω W | ≤ cH n-1 (J u ) n/(n-1) . Since there are finitely many choices of W , |ω| ≤ W |ω W | ≤ cH n-1 (J u ) n/(n-1) . This concludes the proof of (i).

We finally turn to the estimate for the trace. We denote by Γ = {(x ′ , ϕ(x ′ )) : x ′ ∈ (-1/2, 1/2) n-1 } the part of the graph of ϕ which is contained in ∂Ω int , and define for ξ ∈ S ω Γ ξ = {x ∈ Γ :

H 1 (R x ξ ∩ ω * ξ ) > 0} ,
which obeys Π ξ ω Γ ξ ⊂ Π ξ ω * ξ and therefore H n-1 (ω Γ ξ ) ≤ cH n-1 (J u ). We denote by T u : Γ → R n the trace of u on Γ. For any x ∈ Γ \ ω Γ ξ we have Proof of Theorem 1. We choose finitely many sets Ω 1 , . . . , Ω M , Ω int 1 , . . . , Ω int M of the type entering Proposition 6 (up to rotations and translations) such that ∂Ω is covered by the Ω int i . Then we choose finitely many cubes Ω int M +1 , . . . , Ω int K which cover Ω\∪Ω int i such that the double cubes Ω M +1 , . . . , Ω K are contained in Ω. On each of them we can apply Proposition 6 (using ϕ(x ′ ) = r for the cubes). We can assume that J u is sufficiently small that c L H n-1 (J u )

n/(n-1) ≤ 1 3

min{|Ω int i ∩ Ω int j | : Ω int i ∩ Ω int j = ∅} ,
otherwise ω = Ω will do (here c L is the constant in Proposition 6, which depends only on p, n and the Lipschitz constant of Ω). We obtain finitely many exceptional sets ω 1 , . . . , ω K and ω 1 Γ , . . . , ω M Γ (for the inner cubes we do not need the boundary estimate) and finitely many affine functions a 1 , . . . , a K . These functions are in a finite-dimensional space, and since Ω is connected i,j

a i -a j L 1 (Ω int i ∩Ω int j )
is a seminorm on the space of the K-uples of affine functions (a 1 , . . . , a K ), which vanishes on the subspace {a 1 = • • • = a K }. By the assumption on J u also i,j a ia j L 1 (Ω int i ∩Ω int j \(ω i ∪ω j ))

is a seminorm. Therefore, there exists c > 0 such that inf a affine

sup i a i -a L 1 (Ω i ) ≤ c i,j a i -a j L 1 (Ω int i ∩Ω int j \(ω i ∪ω j ))
and we can choose a unique affine function. The constant depends only on Ω because for any α there is c α > 0 such that a L 1 (Ω) ≤ c α a L 1 (E) for all a affine, E ⊂ Ω with |E| > α.

Choosing α as the minimum of the nonzero |Ω int i ∩ Ω int j |/3, then for the same indices

|Ω int i ∩ Ω int j \ (ω i ∪ ω j )| ≥ |Ω int i ∩ Ω int j | -|ω i | -|ω j | ≥ α.

Figure 1 :

 1 Figure1: Sketch of the geometry. The jump set is short with respect to the diameter of the domain Ω and cannot cut it into two "large" blocks. However, small regions can be disconnected from the rest, either inside (a) or close to the boundary (b); also in regions which are "almost" disconnected from the majority region (c) one cannot obtain a uniform estimate.

Figure 3 :

 3 Figure 3: Two possible difficulties in Proposition 2. (a):In choosing the simplex we need to make sure that the jump set does not intersect its edges. (b): Since the jump set will "shadow" points y with [y, z i ] ∩ J u = ∅, we need to make sure that we choose the simplex so that no larger part of J u is close to its vertices. Analogously e(u) should not concentrate on the edges of the simplex and around the vertices.

  ) give Q\ω |w| p (y)dy ≤ cH(z 0 ) ≤ c e(w) p L p (Q 1 ) ≤ c e(u) p L p (Q 1 )

Figure 4 :

 4 Figure 4: Sketch of the geometry in Proposition 3.

Figure 5 : 0 ξ

 50 Figure 5: Sketch of the construction in Proposition 6.

  |ξ • u| p (xtξ) ≤ c R x ξ f p dH 1 for a.e. t ∈ (0, 1/2)which implies, by the property of traces,|ξ • T u| p (x) ≤ c R x ξ f p dH 1 for H n-1 -a.e. x ∈ Γ \ ω Γ ξ .Let now V ⊂ S be a set of n vectors, which are automatically linearly independent, and setω Γ = ∪ ξ∈V ω Γ ξ . Then integrating over x ∈ Γ \ ω Γ gives Γ\ω Γ |T u| p dH n-1 ≤ c V ξ∈V Γ\ω Γ ξ |ξ • T u| p dH n-1 ≤ c e(u) p L p (Ω) ,which is the required estimate.

(i) There is a set ω ⊂ Ω int such that L n (ω) ≤ c L (H n-1 (J u )) n/(n-1) (4.3) and ua L np/(n-1) (Ω int \ω) ≤ c L r 1-1/p e(u) L p (Ω) .

(4.4)

and the trace T u of u on Γ obeys

The constant c L depends only on p, n and the Lipschitz constant L of ϕ.

By affine isometry we mean a map

The geometry is illustrated in Figure 5.

Proof. The first part of the proof is similar to the proof of Proposition 3. Replacing u by ũ(x) = A -1 u(rAx) we can assume r = 1 and A(x) = x. By Proposition 2 applied to the cube Q = (-1, 1) n there are a set ω 0 ⊂ Q and an affine function a with e(a) = 0 such that |ω 0 | ≤ cH n-1 (J u ) and

We can assume without loss of generality a = 0, otherwise we replace u by ua. We let

Then the following holds: (i) the outer normal

(ii) for any ξ ∈ S and any x ∈ Ω int the line x + Rξ intersects the surfaces {(x ′ , x n ) :

For ξ ∈ S and x ∈ Ω int we define the ray R x ξ = (x + Rξ) ∩ Ω , (4.5) see Figure 5. The definition of S ensures that R x ξ is a segment and that H 1 (R x ξ ∩ Q) ≥ 2.


Michael J Norman

C Chandre

T Uzer

Peijie Wang

Nonlinear dynamics of ionization stabilization of atoms in intense laser fields

come

I. INTRODUCTION

"Stabilization" is the term used for the counterintuitive behavior of atoms in which increasing the laser intensity does not lead to increased ionization -on the contrary, it may lead to decreased ionization in some regimes. Stabilization has been the subject of so many publications that even a cursory review is beyond the scope of our manuscript (for some early references, see Refs. [1][2][3][4][5][6][7][8], or Ref. [9] for a review).

The purpose of our manuscript is to provide new insights into the stabilization phenomenon, gained by viewing it through the unconventional point of view of nonlinear dynamics. This is not to say that the use of classical mechanics to examine stabilization is a novelty: A number of stimulating studies were performed during the time frame mentioned above (e.g., Ref. [6]). However, the intervening two decades have seen the realization of lasers with intensities far beyond what was expected at the time, and quite a few phenomena such as recollisions [10][11][12] have been discovered, requiring the inception of various theoretical models and techniques for understanding ionization (or the lack of it) in intense laser fields.

In the intervening decades laser pulses have become ultrashort (down to about 100 attoseconds). The general idea is to go down to the timescale of the electron to capture its dynamics. Such ultrashort timescales were not part of stabilization research at the time, which typically considered fairly long pulses of several tens of laser cycles. Can any insights about such a long-pulse phenomenon be relevant for understanding the motion of the electron on its timescale? As we will show below, no matter how long the pulse is, the fate of the electron, and hence stabilization, is sealed quite early in the pulse, and specifically during a few first few laser cycles. In that sense, stabilization turns out to be a short-timescale phenomenon.

In what follows we will be answering two questions: First, what accounts for the lack of ionization in ultraintense laser fields, and secondly, what causes ionization to decrease with increasing intensity? It turns out that the answers to these questions, and the hidden short-time nature of the phenomenon, are governed by the same un-derlying mechanism, namely the periodic orbits of the atomic electron in the laser field and the phase-space structures associated with these periodic orbits. We will show that electron trajectories which can reach the vicinity of the periodic orbit behave completely differently from those trajectories that fail to do so, thereby accounting for the observations connected with stabilization.

In the process of answering these questions we will also address some highly practical aspects of stabilization, namely :

1) How does one decide whether an atomic electron has ionized in an intense laser field? More specifically, are the so-called distance and energy criteria for ionization interchangeable?

2) What role do the ramp-up, plateau and ramp-down of the laser field play in stabilization?

In Sec. II we start by performing the numerical experiment of Ref. [6] and assess the importance of the choice of the ionization criterion. In Sec. III we provide an explanation for the lack of complete ionization at high intensity and for the decrease of the ionization probability as the laser intensity is increased. This explanation is based on a periodic orbit analysis of the electronic dynamics. We investigate the role of each part of the laser pulse, namely the ramp-up, the plateau and the ramp-down, in the stabilization process.

II. IONIZATION STABILIZATION

We consider the following classical model for a onedimensional single active electron atom interacting with an intense linearly polarized laser field in the dipole approximation :

H(x, p, t) = p 2 2 - 1 √ x 2 + 1 -E 0 f (t)x sin ωt, ( 1
)
where p is the momentum of the electron canonically conjugate to the position x along the polarization axis, E 0 is the amplitude of the electric field, f (t) is the pulse envelope (schematically represented on Fig. 1) and ω is the frequency. We have chosen a soft Coulomb potential [13,14] as the interaction potential between the elec- tron and its parent ion. In Ref. [15], the impact of the choice of potential on stabilization has been analyzed. In all the computations presented in this article, the laser frequency is ω = 0.8 a.u. which corresponds to a wavelength of 57 nm. We begin with revisiting a numerical experiment performed in Ref. [6]. For the envelope (see Fig. 1), we select a smooth turn-on consisting of a linear ramp-up of duration T u = 6T , a plateau of T p = 44T , where T = 2π/ω is the period of the laser, and an abrupt ramp-down, i.e., of duration T d = 0. We consider a large ensemble of trajectories which are initially bounded, i.e., whose energy, defined as the sum of the kinetic energy plus the soft Coulomb potential, is negative. Typically for a given E 0 the number of initial conditions we consider is on the order of several millions in order to have good statistics. As a function of the amplitude E 0 , we measure the ionization probability where ionization is determined using the energy criterion as in Ref. [6]. In Fig. 2, the ionization probability is represented as a function of the amplitude of the electric field E 0 (bold gray curve, red online). Contrary to Ref. [6] we observe a roughly monotonic increase of the ionization probability with E 0 as it could be naively guessed (the more intense the laser field is, the more ionization there is). Now, instead of using the energy criterion as ionization criterion, we use a distance criterion. This distance criterion appears to be more subjective since it depends on an arbitrarily chosen threshold. Here the chosen threshold is 50 a.u.; we have checked that the results are not changed qualitatively by changing this threshold significantly. The resulting ionization curve is depicted in Fig. 2 (thin gray curve, red online). Using this distance criterion we are able to FIG. 2. Ionization probabilities as functions of E0 obtained for Hamiltonian (1) with a ramp-up of duration Tu = 6T , a plateau of duration Tp = 44T . The bold curve (red online) is without a ramp-down (T d = 0) and using the energy criterion for ionization. The thin gray curve (red online) is without a ramp-down (T d = 0), using a distance criterion with a threshold of 50 a.u.. The black curve is with a ramp-down of duration T d = 6T and using the energy criterion.

reproduce, at least qualitatively, the surprising results obtained in Ref. [6], namely the lack of complete ionization at very high intensity and the global decrease of the ionization probability as the intensity of the laser field is increased (see Fig. 2). Here we first notice the importance of the ionization criterion to observe the ionization stabilization phenomenon. We performed a third numerical experiment by including a ramp-down with T d = 6T and by using the energy criterion. Contrary to the case with an abrupt ramp-down (T d = 0), the energy criterion leads to stabilization (see Fig. 2, black curve). From these numerical simulations, it appears that the energy criterion can only be used in the presence of a ramp-down, whereas the distance criterion always displays the stabilization phenomenon. In what follows, we only use the distance criterion unless specified otherwise. Below we provide an explanation for the differences resulting from the two choices of ionization criterion. This is a direct result of the mechanism by which an electron remains trapped near the parent ion at a very high intensity, a mechanism we unravel using nonlinear dynamics.

III. ANALYSIS OF STABILIZATION FROM THE NONLINEAR DYNAMICS PERSPECTIVE A. Role of the plateau

We first elucidate the role of the plateau by providing the scenario for the lack of complete ionization at high intensity. Since the amplitude of the electric field is large, we consider first the Strong Field Approximation (SFA) which consists of neglecting the soft Coulomb interaction between the electron and the parent ion. In SFA, Hamiltonian (1) during the plateau becomes

H SFA (x, p, t) = p 2 2 -E 0 x sin(ωt + φ), (2)
where we have included a laser phase to take into account the ramp-up, i.e., we consider that the plateau starts at t = 0. This Hamiltonian is integrable and the trajectories can be explicitly provided :

p(t) = p 0 - E 0 ω (cos(ωt + φ) -cos φ) , x(t) = x 0 + p 0 + E 0 ω cos φ t - E 0 ω 2 (sin(ωt + φ) -sin φ) ,
where x 0 and p 0 are the initial position and initial momentum (at t = 0, i.e., at the beginning of the plateau). We notice a linear drift with velocity v d = p 0 + (E 0 /ω) cos φ. If this drift velocity is non-zero, the electron will eventually move far away from the parent ion, provided that the plateau is sufficiently long. We notice a very special set of initial conditions for which the drift velocity vanishes. This set in phase space is composed of an infinite family of (parabolic) periodic orbits with the same period as the laser field, related to each other by translation along the x-axis due to the continuous symmetry of Hamiltonian (2). More specifically, these periodic orbits have the equations

p(t) = - E 0 ω cos(ωt + φ), (3a)
x(t) = x 1 - E 0 ω 2 sin(ωt + φ), (3b)
where x 1 is an arbitrary position, the position of the center of the periodic orbit. If the electron is on one of these orbits at the beginning of the plateau, it will stay on it for an infinite amount of time. These invariant structures are the keystone to understand the lack of complete ionization at high intensity.

Even though this set of periodic orbits is of measure zero and therefore it could be argued that it does not contribute significantly to the ionization probability, the nonlinearity introduced by the soft Coulomb interaction provides the reasons for the importance of this set of periodic orbits in the stabilization process, as we show next.

The question is what happens to the rather simple SFA picture of the dynamics when the Coulomb field is taken into account? As the effect of the soft Coulomb potential increases, e.g., by increasing an effective charge up to 1, most of the periodic orbits (3) are broken, except a finite number of them. Three periodic orbits are of particular interest and organize the dynamics. All three periodic orbits have the same period, the period of the laser field. One of them is symmetric with respect to the x = 0 axis, centered at (x, p) = (0, 0) and weakly hyperbolic. This periodic orbit is denoted by O. The other two periodic orbits of interest are elliptic and symmetric with each other, centered approximately at (x, p) = (±E 0 /ω 2 , 0). These orbits are denoted O ± . The property of a periodic orbit to be hyperbolic, parabolic or elliptic features refers to the linear stability analysis (see Ref. [START_REF] Cvitanović | Chaos: Classical and Quantum[END_REF] for more detail). In general, an elliptic periodic orbit is stable and the motion around it is similar to the motion around a stable equilibrium point of a simple pendulum, whereas an hyperbolic orbit is (linearly) unstable, in the sense that almost all small perturbation around it will drive the motion away from it (along its unstable manifold). The chaotic behavior of the system originates in the neighborhood of hyperbolic periodic orbits.

Figure 3 displays the three periodic orbits of Hamiltonian (1) with f (t) = 1 and for E 0 = 5. At the high values of the intensity we consider in this paper, O is very close to the periodic orbit given by Eqs. (3) for x 1 = 0 (see Fig. 3 where it is difficult to distinguish the orbit (3) from the orbit O). In Fig. 4 we represent the maximum eigenvalue of the monodromy matrix associated with O as a function of E 0 (see Ref. [START_REF] Cvitanović | Chaos: Classical and Quantum[END_REF] for more details). We notice that for E 0 small, O is elliptic and, as E 0 is increased, undergoes a bifurcation around E 0 = 1.34. The orbit O turns hyperbolic and two asymmetric and elliptic periodic orbits are created, these being O ± . The basic structure of phase space is a horizontal 8-shape with two symmetric elliptic islands with the orbits O ± in their centers and the orbit O at the crossing. A Poincaré section of selected trajectories consisting of a stroboscopic plot of period T is represented on Fig. 5 for E 0 = 1.5, i.e., close to the bifurcation. This picture persists at higher value of E 0 , although complicated by the creation of a dense chaotic tangle around these orbits. A Poincaré section of the region close to the weakly periodic orbit is depicted in Fig. 6 for E 0 = 5. The dense chaotic region is clearly observed, extending to high values of the position of the electron. Trajectories starting in the white region ionize very quickly, contrary to the trajectories started in the blue region which can be potentially trapped for an arbitrary long time. The two symmetric elliptic islands are present at this very high value of the intensity (see Fig. 7), demonstrating the stabilization effect caused by the absence of drift velocity in the SFA approximation. In this chaotic region, there are many other periodic orbits and some of them might even be elliptic (depending on E 0 ) forming islands of stability (trapping trajectories at all times). It should be noted that the two periodic orbits O ± remains elliptic for all values of E 0 , becoming closer and closer to parabolic as E 0 is increased. In Fig. 8, we represent Greene's residue [START_REF] Greene | [END_REF]18] of O ± as a function of E 0 . We recall that that if the residue is between 0 and 1, the periodic orbit is elliptic. This property is very interesting since the excited states around O ± (which extends up to the quiver radius, see Fig. 3) are stable states.

In a nutshell, the effect of the soft Coulomb interaction is to create a chaotic layer which strongly affects the transport properties of the system, trapping trajectories for a sufficiently long time compared to the duration of the plateau. These trajectories will not ionize, in the sense that they will not significantly depart from a region close to the core, contrary to the trajectories experiencing a drift velocity. This chaotic region explains the lack of complete ionization, even at very high values of the laser intensity. Now, concerning the results obtained using the energy criterion, we notice that by looking at the energy of the non-ionizing trajectories (where nonionizing feature has been determined using a distance criterion), given the of the momentum which extend up to E 0 /ω, most of these trajectories have positive energy. For instance, we have displayed in red on Fig. 3 the points of the orbits O and O ± which have negative energy (and would be considered as non-ionizing using an energy criterion). Therefore, according to an energy criterion, most of the trapped trajectories are considered as ionizing even though they are of a very different nature than the ones experiencing a drift velocity. Globally the effect of the plateau is to discriminate between the trajectories with a drift velocity with the ones trapped in the chaotic tangle. The longer the pulse is, the more trapped trajectories leave the chaotic tangle, and eventually these will be considered as ionized (even using a distance criterion). In Fig. 9, we represent the ionization probabilities as functions of E 0 for different values of the plateau duration. It confirms that the ionization probability increases with the plateau However it should be noticed that due to the presence of elliptic islands, the ionization will never be complete for most values of the laser intensity, since the trajectories in these regions are trapped for all times, preventing the ionization probability to reach 1 regardless of the length of the plateau. It should be noticed that the finer structure of the ionization curve, like the oscillations, are qualitatively similar (located at approximately the same values of E 0 ), regardless of the duration of the plateau.

Next we look at the size of the chaotic/trapping region as a function of E 0 . In order to do this, we launch a high number of trajectories in a window around the periodic orbit in the Poincaré section. The window we consider is

(x, p) ∈ [-20, 20]×[-E 0 /ω-3/2, -E 0 /ω+3/2],
where we recall that the periodic orbit O is locate approximately at (0, -E 0 /ω). We look at trajectories which have left this window after 20 and 50 laser cycles. In Fig. 10, we represent the probability of a trajectory to leave the window of interest. We compare this probability with the ionization probability of Fig. 2 with a ramp-up of duration T u = 6T and no ramp-down. The probability to leave the window exhibits the main features of the ionization probability. First it is decreasing (meaning that the size of the chaotic region increases with E 0 ), therefore displaying the stabilization process. In addition, we notice that this curve displays the same finer structures (the oscillations as function of E 0 ). Therefore the structure of the ionization curves of Fig. 2 is a nonlinear effect caused by the presence of both the soft Coulomb potential and the laser field. As expected the ionization probability is smaller due to the fact that the average intensity felt by the electron is smaller when there is a ramp-up. The ramp-up has also another impact on the ionization probability and this will be analyzed below. The global decrease of the ionization curves is due to the fact that, as E 0 is increased the system becomes closer to integrable as shown in Figs. 4 and8 where the periodic orbits becomes closer to parabolic.

It should be noted that the same periodic orbits O and O ± which control the stabilization process, also control the recollision process by allowing an electron far away from the parent ion to come back to it and recollide. These periodic orbits were coined recolliding periodic orbits (RPO) in Refs. 20].

B. Role of the ramp-up

The role of the ramp-up is important in setting up the right conditions for some of the trajectories to be put inside the chaotic/trapping region at the beginning of the plateau. In order to illustrate this argument, we first consider the SFA approximation, by considering the following Hamiltonian :

H u (x, p, t) = p 2 2 -E 0 x t T u sin(ωt + φ),
where here the phase φ is used to take into account the transient time when the SFA approximation is not valid. The momentum and the position at the end of the rampup are

p u = p 0 - E 0 ω cos(ωT u + φ) - sin(ωT u + φ) -sin φ ωT u , x u = x 0 - E 0 ω 2 (sin(ωT u + φ) + sin φ +2 cos(ωT u + φ) -cos φ ωT u .
If the duration of the ramp-down is an integer multiple of the laser period then it reduces to which is a point on a periodic orbit given by Eqs. ( 3) with x 1 = x 0 -(E 0 /ω 2 ) sin φ, provided that p 0 = 0. If the initial distribution at the beginning of the ramp-up is close to zero then the role of the ramp-up (at least in the SFA approximation) is to promote these initial conditions to a region where periodic orbits organize the dynamics and prevent ionization. Taking into account the soft Coulomb potential, the role of the ramp-up is to promote some of the trajectories to the chaotic/trapping region organized by the three periodic orbits O and O ± .

p u = p 0 - E 0 ω cos φ, x u = x 0 -2 E 0 ω 2 sin φ,
In order to demonstrate this effect of the ramp-up, we depict the position of an ensemble of electrons in phase space at different stages in the ramp-up. Initially these trajectories are bounded, in the sense that their energy (sum of kinetic energy and soft Coulomb potential is negative). In Fig. 12, four stages are represented : at t = 0, t = 2T , t = 4T and t = 6T . We notice that initially none of the trajectories are on the periodic orbit. However, as early as half of the ramp-up, some of the trajectories are already on the periodic orbit or nearby. If the ramp-up is not sufficient to promote the trajectories to the chaotic/trapping region, for instance, if it is too short (shorter than half a laser cycle for instance), then all the trajectories will experience a significant drift velocity and will be all ionized. Therefore the ramp-up does not play a role in the stabilization mechanism but allows part of the initial conditions to undergo this stabilization.

C. Role of the ramp-down

Using the energy criterion, the presence of the rampdown seems essential for the stabilization. However we have shown above that the mechanism for stabilization has nothing to do with the ramp-down. In this section we analyze the effect of the ramp-down, showing that it only plays a role in the dynamics. However it is essential if one wants to use the energy criterion for ionization. In order to illustrate this, we consider the effect of the ramp-down on the prototypical trapped trajectory in the SFA approximation, namely the periodic orbit given by Eqs. (3). We shift the origin of time, so that the rampdown starts at t = 0 (and we assume that the ramp-up and plateau are integer multiples of the laser period). The duration of the ramp-down is T d . We consider the following Hamiltonian in the SFA approximation during the linear ramp-down:

H d (x, p, t) = p 2 2 -E 0 x T d -t T d sin(ωt + φ).
At t = 0 (i.e., at the end of the plateau), we assume that the trajectory is on the periodic orbit, i.e., the initial condition is x 0 = -(E 0 /ω 2 ) sin φ and p 0 = -(E 0 /ω) cos φ at the beginning of the ramp-down. At the end of the ramp-down, the final momentum is

p f = E 0 ω sin φ -sin(ωT d + φ) ωT d .
If the duration of the ramp-up is an integer multiple of T , then the final momentum vanishes. Regardless of the distance of the electron to the parent ion, its energy is then negative. Therefore the electron located on the SFA periodic orbit at the end of the plateau is considered ionized at the end of the ramp-down according to the energy criterion. However if the duration of the ramp-down is not an integer multiple of T , the final momentum can be large (up to E 0 /ω) and consequently can be considered as ionizing according to the energy criterion. It should be noted that using the distance criterion, the effect of the ramp-down can be neglected; all the ionization probability curves are nearly identical. In Fig. 13 we represent the ionization probabilities obtained with an energy criterion for different durations of the ramp-down, one of them being not an integer multiple of the laser period. Provided that the duration of the ramp-up is an integer multiple of the laser period, the two curves at T d = T and T d = 3T demonstrate little effect of the ramp-down, even using the energy criterion as ionizing criterion.

IV. CONCLUSION

We have shown that the main mechanism for stabilization of atoms in a strong laser field is driven by a set of periodic orbits. The counter-intuitive nature of this phenomenon can be explained by the fact that the orbits mostly lay outside the bounded region. We have elucidated the role of each phase in the laser pulse in the stabilization process. During the turn-on of the laser field, when the Coulomb and laser fields are of comparable strength, an electron is pushed towards a dense chaotic tangle in the neighborhood of these periodic orbits, trapping some trajectories for an arbitrarily long time (even for infinite time depending on the initial conditions). The role of the plateau is to seal the fate of the electron by discriminating the ones that are non-ionizing from the ionizing trajectories. We have showed that the role of the ramp-down is rather minor, allowing the use of the energy criterion, whereas, in the absence a ramp-down a distance criterion has to be used to detect the stabilization phenomenon.

FIG. 1 .

 1 FIG.1. Schematic representation of the pulse envelope f (t) (left axis, bold line) and electric field E(t) = -E0f (t) sin ωt (right axis, thin line) used in Hamiltonian(1). Since the plateau is much longer that the ramp-up and the ramp-down, it has not been fully depicted.

FIG. 3 .

 3 FIG. 3. Projection of the three periodic orbits O and O± of Hamiltonian (1) with f (t) = 1 and E0 = 5. The light gray orbit is O+, the darker gray orbit is O-and O is represented in black. The gray parts (red online) of the orbits are where the energy (defined as the sum of the kinetic energy and soft Coulomb potential) is negative. The thin line corresponds to the periodic orbit in the SFA approximation given by Eq. (3), almost indistinguishable from O.

FIG. 5 .

 5 FIG. 5. Poincaré section (stroboscopic plot with period T ) of trajectories of Hamiltonian (1) for E0 = 1.5. The positions of the periodic orbits are indicated by a cross for O and by circles for the elliptic periodic orbits O±.

FIG. 7 .

 7 FIG. 7. Inset of Fig. 6.

FIG. 9 .

 9 FIG.9. Ionization probabilities as functions of E0 for different durations of the plateau. From black to lighter gray (red online), the total duration of the pulse is 50, 100 and 200 laser cycles, with a 6 laser cycle ramp-up and no ramp-down.

FIG. 10 .

 10 FIG. 10. Probability of a trajectory of Hamiltonian (1) with no ramp-up and no ramp-down to leave the window (x, p) ∈ [-20, 20]×[-E0/ω-3/2, -E0/ω+3/2] after 20 (light gray, red online) and 50 (dark gray) laser cycles. For comparison, we depict the ionization curve of Fig. 2 obtained with a ramp-up of duration Tu = 6T and no ramp-down (black curve).

FIG. 11 .

 11 FIG.11. Ionization probabilities as functions of E0 for different durations of the ramp-up. From black to lighter gray (red online), the duration of the ramp-up is T /10, T /2, T and 6T with a 44 laser cycle plateau and no ramp-down.

FIG. 12 .

 12 FIG. 12. Positions of an ensemble of trajectories of Hamiltonian (1) for E0 = 5 at different stages of the ramp-up, t = 0 (upper left panel), t = 2T (upper right panel), t = 4T (lower left panel), and at the end of the ramp-up t = 6T (lower right panel). The periodic orbit O is depicted in gray (red online).

FIG. 13 .

 13 FIG. 13. Ionization probabilities as functions of E0 obtained from Hamiltonian (1) for different durations of the rampdown. From black to lighter gray (red online), the duration of the ramp-down is T d = T /5 (black curve), T d = T (dark gray, maroon online) and T d = 3T (light gray, red online), with a 6 laser cycle ramp-up and a 44 laser cycle plateau. The energy criterion is used to determine ionization.

ACKNOWLEDGMENTS

T.U. acknowledges R. Grobe for some clarifications. M.J.N and T.U. acknowledge funding from the NSF. The research leading to these results has received funding from the People Program (Marie Curie Actions) of the European Unions Seventh Framework Program No. FP7/2007-2013/ under REA Grant No. 294974. P.W. acknowledges support from the Scientific Research Base Development Program of the Beijing Municipal Colleges Higher Education Talent Programs (No. 067135300100).


F Grognard

Y Chitour
email: yacine.chitour@math.u-psud.fr

G Bastin
email: bastin@auto.ucl.ac.be

EQUILIBRIA AND STABILITY ANALYSIS OF A BRANCHED METABOLIC NETWORK WITH FEEDBACK INHIBITION

Keywords: Metabolic network, small gain, stability, limit cycle

This paper deals with the analysis of a metabolic network with feedback inhibition. The considered system is an acyclic network of mono-molecular enzymatic reactions in which metabolites can act as feedback regulators on "preceding" enzymes of their own pathway, and in which one metabolite is the root of the whole network. We show, under mild assumptions, the uniqueness of the equilibrium. In the simplified case where inhibition only acts on the reactions having the root of the network as substrate, we show that the equilibrium is globally attractive. This requires that we impose conditions on the kinetic parameters of the metabolic reactions.

INTRODUCTION

The cellular metabolism is defined as the (huge) set of biochemical reactions that occur inside a living cell for growth and reproduction. It is usually represented by an intricate network connecting the involved biochemical species (called "metabolites"). The pathways of the network are called "metabolic pathways". In the metabolic engineering literature, it is widely accepted that "despite their immense complexity, metabolic systems are characterized by their ability to reach stable steady states" (quoted from [START_REF] Stephanopoulos | Metabolic Engineering: Principles and Methodologies[END_REF], Chapter 4). It should however be fair to recognize that a mathematical analysis of this fundamental stability property is a difficult question which was not much investigated. We shall limit ourselves to simple metabolic pathways which are made up of sequences of mono-molecular enzyme-catalysed reactions in the form X s → X p .

Those reactions can be inhibited by the presence of other metabolites in the network. Without this inhibition, the stability analysis is straightforward; we will then concentrate on networks with inhibition, like the one given by the aspartate amino-acid pathways [START_REF] Umbarger | Amino acid biosynthesis and regulation[END_REF], see Figure 1. In this network, each produced amino-acid inhibits an enzyme of its own pathway. This action can be seen as a negative feedback, that regulates the behavior of the network. Indeed, if we, for example, consider a large excess of isoleucine (X 20 ), the reaction X 16 -→ X 17 is shut down, so that the concentration of isoleucine is progressively reduced.

In Section 2, a model of metabolic networks such as the one of Figure 1 will be presented. The equilibria of these models will then be studied in Section 3, followed by a stability analysis in Section 4, where global attractivity of a single equilibrium is shown. x 20 Fig. 1. Metabolic network representing the aspartate amino-acid pathways: the solid lines represent the reactions and the dash-dotted lines the inhibition produced by the state at the start of the arrow onto the reaction that lies at the end of the arrow. The root of the metabolic pathway is x 1 (aspartate), and the products are the corresponding amino acids: lysine (x 9 ), methionine (x 14 ), threonine(x 16 ), and isoleucine (x 20 ).

The non-genericity of the stability of the equilibrium is then illustrated in Section 5.

MODEL OF A METABOLIC NETWORK

In our model of a metabolic network made of enzymecatalysed reactions in the form X s → X p which are inhibited by other metabolites of the network, we will denote the inhibition factor by x [p] , a n p dimensional vector containing the molar fractions of all the metabolites inhibiting the reaction X s → X p . This reaction is characterized by a velocity ϕ sp (x s , x [p] ).

On the other hand, some metabolites of the considered network are used in stages of the metabolism that are not modeled in the considered network. Therefore, we must include consumption terms in the model for those reactions in the form X s → ... We will generically denote those terms ϕ s0 (x s ).

We impose that these reaction velocities satisfy the following assumption:

Assumption 1.

• For all s, p such that the reaction X s → X p belongs to the metabolic network, the function ϕ sp (x s , x [p] ) satisfies ϕ sp (0) = 0, is non-decreasing in x s for x s ≥ 0 and non-increasing in x

[p] j . Note that this must also be valid when the value of p is 0.

In order to define the class of metabolic networks that we consider, we need the following definition from graph theory: Definition 1. A directed graph is called an arborescence if, from a given node x, known as the root node, there is exactly one elementary path from this node to any other node y.

The metabolic networks that we consider in this paper then satisfy the following assumption Assumption 2.

• the involved species are denoted

X 1 , X 2 , • • • , X n
• the graphic representation of the network (with the different metabolites as nodes and the different reactions as oriented edges) is an arborescence with X 1 as root • the inhibition acting on a reaction X s → X p only results from the action of metabolites from the (sub)-arborescence rooted in X p

Stemming from the definition and known properties of arborescences, Assumption 2 has the following consequence on the class of metabolic network that we consider:

(i) Each metabolite is produced by a single other metabolite; (ii) There is no cycle of reactions; With these definitions and notations, we shall now define a mass-balance dynamical model in the form

ẋ = Φ(x) -µx + ce 1 where x = (x 1 , • • • , x n ) T ∈ IR n ,
and x i denotes the molar fraction of the metabolite X i inside the cell. The factor µ ≥ 0 represents the specific growth rate of the cell: we assume that the cell metabolism is analyzed during a period of exponential cell growth with a constant specific growth rate µ. The vector e 1 = (1, 0, • • • , 0) T and the scalar c denote the constant supply rate of the metabolite X 1 at the root of the network. The function Φ includes all the reaction velocities, whether they correspond to reactions inside the network or reactions consuming the metabolites of the network for use in subsequent stages of the metabolism.

In order to specify Φ(x), we introduce the following notations:

Notation 1.

• P(j) = {k| the reaction X j → X k belongs to the network }. P(j) defines the set of all metabolites that are produced by reactions having X j as substrate. If there is a consumption term in the form ϕ j0 (x j ) in the derivative of x j , the index 0 is included in P(j).

• A(j) = {k|X k belongs to the arborescence with its root in X j }. "0" is not included in A(j).

It can easily be seen that P(j) \ {0} is a subsets of A(j) under Assumption 2.

From the arborescence structure, it is clear that we can separate the metabolites into three different families:

• the root X 1 : we suppose that there is a constant supply rate c of X 1 so that the corresponding mass-balance equation is the following:

ẋ1 = c - k ∈P(1) ϕ 1k (x 1 , x [k] ) -µx 1 (1)
• the intermediate metabolite X j , which is the result of the reaction X i → X j :

ẋj = ϕ ij (x i , x [j] ) - k ∈P(j) ϕ jk (x j , x [k] ) -µx j (2)
• the boundary metabolite X j (such that P(j) = {0}), which is the product of a reaction X i → X j :

ẋj = ϕ ij (x i , x [j] ) -ϕ j0 (x j ) -µx j (3)
Under Assumption 2, we can only have

x [j] = x j or x [j] = ∅.
A particular case of this network is the metabolic chain with sequential feedback inhibition (cf. Chitour, [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF])

X 1 → X 2 → • • • → X n .
with the last metabolite X n acting as an inhibitor of the first reaction X 1 -→ X 2 . In this case, all x [j] are empty, except x [2] = x n .

EQUILIBRIUM OF A METABOLIC NETWORK

Based on equations ( 1)-( 2)-(3), we can now compute the mass-balance of the whole arborescence:

d dt ( n l=1 x l ) = c -µ n l=1 x l - {k|0 ∈ P(k)} ϕ k0 (x k ) (4)
and of the arborescence that has its root in

X j d dt ( l ∈A(j) x l ) = ϕ ij (x i , x [j] ) -µ l ∈A(j) x l - {k|0 ∈ P(k) and k ∈A(j)} ϕ k0 (x k ) (5)
Those expressions will be critical in the proof of the following proposition:

Proposition 1. If Assumptions 1 and 2 are satisfied, then:

(A) the system (1)-( 2)-( 3) is positive;

(B) if all ϕ s,p (x s , x [p] ) (where p can be 0) are increasing in x s then there is at most one equilibrium

x = (x 1 , • • • , xn ) of (1)-(2)-(3) in IR n + ; (C) if µ > 0, then system (1)-(2)-(3) has a unique equilibrium x = (x 1 , • • • , xn ) in IR n + .
Moreover, the solutions of ( 1)-( 2)-( 3) are bounded for any initial condition in IR n + .

Proof: (A) is easily seen by considering the system on the boundaries of the positive orthant.

The proofs of (B) and (C) are very similar. We write the proof for (B) and highlight the differences that arise for the proof of (C) .

We will first consider system (2)-( 3) with x 1 = x1 as constant input. For any value of x1 , we will denote by (x 2 , • • • , xn ) the equilibrium of ( 2)-( 3); this equilibrium is a function of x1 , so that we will state that it is

(x 2 , • • • , xn )(x 1
). We will now show, by induction, that every element xi is an increasing function of x1 (resp. non-decreasing in case (C)).

The initial step of the proof considers the equilibrium of (3) with x i = xi as constant input

ϕ ij (x i , x[j] ) -ϕ j0 (x j ) -µx j = 0
where x[j] = xj or x [j] = ∅. When xi = 0, xj = 0 is the only solution. Also, the left-hand side of this equation is an increasing function of xi (resp. nondecreasing in case (C)) and a decreasing function of xj . It is then easily seen that, if we increase xi , xj needs also to be increased (resp.increased or kept constant) to keep this equality satisfied. We then have that, in this case, xj (x i ) is an increasing function such that xj (0) = 0 (resp. non-decreasing function such that xj (0) = 0). When µ = 0 (which can only happen in case (B)), the definition of xj could be limited to an interval [0, xm i ). Let us now make the following induction hypothesis: for a given j, the functions xk (x j ) are increasing (resp. non-decreasing) functions for all k ∈ A(j) with xk (0) = 0. We then study the equilibrium of the massbalance of the arborescence that has its root in X j . From (5):

ϕ ij (x i , x[j] (x j )) -µ l ∈A(j) xl (x j ) - {k|0 ∈ P(k) and k ∈A(j)} ϕ k0 (x k (x j )) = 0
With a similar argument to that of the initial step, we see that xj is an increasing (resp. non-decreasing) function of xi and that xj (0) = 0. The same can be said for all xk with k ∈ A j because they are already increasing (resp. non-decreasing) functions of x j .

By induction, we then have that every xk (x 1 ) is an increasing function of x1 defined on the interval [0, xm 1 ) (resp. non decreasing function defined for all x1 ). An equilibrium of the whole system then has to satisfy the equilibrium of the total mass-balance. From (4), this comes to:

µ n l=1 xl (x 1 ) + {k|0 ∈ P(k)} ϕ k0 (x k (x 1 )) = c
The system admits as many equilibria as this equation has roots. In case (B), the left-hand side is increasing from 0 when x1 increases from 0 to xm 1 (because of the second term). Therefore, if there exists an equilibrium, it is unique. In case (C), the left-hand side is strictly increasing from 0 to +∞ when x1 increases from 0 to +∞ (because of the µx 1 term), so that the equilibrium exists and it is unique.

The final point of (C) is a direct consequence of (4); this implies

d dt ( n l=1 x l ) ≤ c -µ n l=1
x l which clearly implies boundedness of the solutions when µ > 0. ✷ Uniqueness of the equilibrium, especially when it is coupled with boundedness of solutions, gives hope of the possibility of having some general result about the structural global asymptotic stability of the equilibrium. In the next section, we study a special case of feedback inhibition, in which we will be able to prove global attractivity.

STABILITY OF A CLASS OF METABOLIC NETWORKS

In this section, we will study the stability of the equilibrium of a class of metabolic networks belonging to the family that was described in Section 2. In order to do that, we will first present a technical lemma.

Technical lemma

In order to analyze the stability of the considered metabolic networks, we will use the small-gain theorem of De Leenheer, Angeli & Sontag (2003) for the interconnection of monotone systems. Due to space limitation, the details of this theory are omitted here, and we only give two lemmas that we will use in the following sections (and which are proven in [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF]. Let us consider the stability of a system

ẋ = f (x) x ∈ IR n (6)
which is not cooperative, but which satisfies the following assumption:

Assumption 3. Each off-diagonal element of the Jacobian J of f (x) is sign-definite (independent of x).

We then build a new system ẋ = F (x, -v) according to the following construction: for all i, j with j = i, if J ij < 0, then replace x j in f i (x) with the constant -v j . Define u 1 as the vector that contains all the constants v j that are necessary for the preceding construction (not necessarily all j have been required); we denote those constants by v s1 , • • • , v s l . From this construction, it directly appears that the system

ẋ = F (x, -u 1 ) (7)
is cooperative (in the sense defined in De Leenheer et al. (2003)) We then define the output of ( 7) as

y 1 = (x s1 , • • • , x s l )
T . The interconnection of system ( 7) with the trivial system

y 2 = u 2 (8)
through the feedback connections u 1 = -y 2 and u 2 = y 1 results in the original system (6). We can then show the following lemma Lemma 1. Suppose that the solutions of ( 6) are bounded and that system ( 7) has a unique globally attractive equilibrium for any constant input u 1 (this equilibrium defines an input output characteristic ȳ1 = k y1 (u 1 ) for system [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF]). Suppose that the expression of y 1 = k y1 (u 1 ) is unknown, but that it is known that it satisfies equation

y 1 = G(y 1 , u 1 )
then system (6) has a unique globally attractive equilibrium if there exists a norm . such that

sup u1, y1 ∈ I R m ∂G ∂y 1 + ∂G ∂u 1 < 1 (9)
and if the unique equilibrium of system [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF] with u 1 constant is globally attractive for any u 1 .

Proof: see [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF] 4

.2 Feedback inhibition

The metabolic networks that we will consider for our stability study are made up of an arborescence of mono-molecular enzyme-catalysed reactions which satisfy the assumption that only one type of inhibition is present: the inhibition of the reactions that use X 1 as substrate. Also, a metabolite X m can only inhibit the reaction X 1 -→ X p , which is the first reaction of the unique elementary path linking X 1 to X m (this is a consequence of the third point of Assumption 2). This translates into Assumption 4.

• For all k ∈ P(1),

x [k] = ∅.
For simplification of notations, we will denote the elements of P(1)

as {k 1 , • • • , k r }.
A mass-balance model for such a network fits in the structure that was described in the previous section, so that we know that there is a single equilibrium. The model can be written as ( 1)-( 2)-( 3), but is particularized due to the presence of Assumption 4:

ẋ1 = c - r i=1 ϕ 1ki (x 1 , x [ki] ) -µx 1 (10) ẋki = ϕ 1ki (x 1 , x [ki] ) - j ∈P(ki) ϕ kij (x ki ) -µx ki (11) ẋk = ϕ lk (x l ) - j ∈P(k) ϕ kj (x k ) -µx k (12)
where x ki is a product of x 1 and x k of x l = x 1 . We impose the boundedness of the partial derivatives of ϕ ij in the following assumption:

Assumption 5. There exist d ij ≥ 0, α [kj ] b ≥ 0 such that 0 ≤ ∂ϕ ij ∂x i ≤ d ij for all i, j -α [kj ] b ≤ ∂ϕ 1kj ∂x [kj ] b ≤ 0 for all j ≤ r, b ≤ n j
We now define a new notation Notation 2. From the arborescence structure, we know that there exists a unique path from X 1 to any metabolite X s . This path takes the form

X 1 → X kj → • • • X k → X l → • • • → X w → X s ; if X s is
an arbitrary metabolite, we will store the indices of this path (without 1 and s) in C s ; alternatively, if x s corresponds to some x

[kj ] b , we will also denote this path C

[kj ] b . Similarly, we denote g s (k) or g [kj ] b (k) the index of the metabolite that follows X k in the path that connects X 1 to X s . This allows for the following theorem Theorem 2. If Assumptions 1, 2, 4 and 5 are satisfied and

r j=1 n k j b=1    d 1,kj µ + 1 k ∈ C [k j ] b d kg [k j ] b (k) µ + d kg [k j ] b (k)    < µ max ks,c α [ks] c (13)
and ϕ 1ki is bounded for all i ∈ {1, • • • , r} (0 ≤ ϕ 1ki ≤ B ki ), then the equilibrium of system ( 10)-( 11)-( 12) is globally attractive in the positive orthant.

Proof: The proof of this Theorem is an application of Lemma 1 (see [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF]). The decomposition of ( 10)-( 11)-( 12) as in ( 7) is done by replacing all the inhibiting terms x [p] with the constant vector -u 1 .

We see that, despite the fact that the inhibition is classically presented as a negative feedback that regulates the system, we have only been able to prove global attractivity under the restrictive condition (13), while stability is easily seen in the absence of inhibition. This condition is very strong, especially if the specific growth rate is small. We will analyze this condition further in the next section.

LIMIT CYCLES IN METABOLIC NETWORKS

Having obtained the sufficient result for global attractivity of Theorem 2, it is relevant to ask two questions: is condition (13) necessary and sufficient, on the one hand, and is Theorem 2 still valid without condition (13) on the other hand? The answer to the first question is easily seen to be "no": the stability of the equilibrium is retained even if condition ( 13) is slightly violated (a few simulations of simple systems is already convincing). In this section, we will show that the answer to the second question is also "no": without condition (13), the stability can be lost, so that we see that the stability of the metabolic networks is not a simple consequence of the structure of the models. Indeed, in this section, we shall exhibit an example where the equilibrium becomes unstable with a limit cycle (Hopf bifurcation) when condition ( 13) is not satisfied. We will concentrate on the stability of the equilibrium of a simple sequential pathway of four metabolites without branching and with sequential feedback inhibition (that was presented at the end of Section 2). Each metabolite produces a single other metabolite, and X 1 -→ X 2 is inhibited by the last metabolite, X 4 . We can directly apply Theorem 2 to this system:

                   ẋ1 = 1.71 - 1 1 + (x 4 /19) p 3.2x 1 1 + x 1 -0.01x 1 ẋ2 = 1 1 + (x 4 /19) p 3.2x 1 1 + x 1 - 1.4x 2 1 + x 2 -0.01x 2 ẋ3 = 1.4x 2 1 + x 2 - 1.2x 3 1 + x 3 -0.01x 3 ẋ4 = 1.2x 3 1 + x 3 - x 4 1 + x 4 -0.01x 4 (14)
where we take ϕ 1 (x 1 , x 4 ) = which is Hurwitz for p < 56.4519, and is not Hurwitz for p larger than that value. This transition from a stable to an unstable equilibrium is illustrated on Figure 2, where the time responses of the four states is illustrated for the value of p = 0 (no inhibition), p = 10 (weak inhibition), and p = 60 (strong inhibition). In the latter case, oscillations appear. This corresponds to a limit cycle in the state-space. A Hopf bifurcation has taken place in p = 56.4519. Despite this oscillation, the reaction rates for the three reactions X 2 -→ X 3 , X 3 -→ X 4 , and X 4 -→ ... are close to their maximum after the transient. The only limiting reaction is X 1 -→ X 2 which, due to the inhibition is far from its maximum reaction rate.

CONCLUSION

In this paper, we have shown that a large class of models of metabolic system only has a single equilibrium. We then have proven that, under a small gain condition, this equilibrium is globally attractive in the particular case where inhibition only acts on reaction having the root as substrate. Finally, we have shown that stability of this equilibrium is not a generic property of the metabolic systems: a condition needs to be imposed on the parameters to have stability (similar to the small gain condition that we found).

Fig. 2 .

 2 Fig.2. Evolution of the states of system (14) for p = 0 (dotted line), p=10 (dash-dotted line) and p = 60 (solid line).

  At the equilibrium, we can see that the characteristics polynomial of the Jacobian matrix has the form (s + 0.0140)(s + 0.0135)(s + 0.0130)(s + 0.0125) +4.2 10 -7 p(s + 0.01)

					1 1+(x4/19) p	3.2x1 1+x1 . The
	parameter p mainly influences the maximal slope of
	the inhibiting factor	1 1+(x4/19) p (which takes place in
	x 4 = 19). Condition (13) becomes
	α µ	d 1 µ	+ 1	d 2 µ + d 2	d 3 µ + d 3

< 1 ≡ p < 3.16 10 -5 which is a condition very similar to what was obtained in Chitour et al. (2003). This condition is very strong because it is based on a small-gain analysis. System (14) has a single equilibrium in x = (19, 19, 19, 19) T .