Michele Correggi

Nicolas Rougerie
email: nicolas.rougerie@lpmmc.cnrs.fr

EFFECTS OF BOUNDARY CURVATURE ON SURFACE SUPERCONDUCTIVITY

We investigate, within 2D Ginzburg-Landau theory, the ground state of a type-II superconducting cylinder in a parallel magnetic field varying between the second and third critical values. In this regime, superconductivity is restricted to a thin shell along the boundary of the sample and is to leading order constant in the direction tangential to the boundary. We exhibit a correction to this effect, showing that the curvature of the sample affects the distribution of superconductivity. Contents 1. Introduction and Main Result 1 2. Correction to the 1D Energy 6 3. Estimates of the Energy Density 10 4. From the Energy Density to the Order Parameter 13 References 17

Introduction and Main Result

The response of type-II superconductors to external magnetic fields is a rich source of fascinating mathematical problems [START_REF] Béthuel | Progress in Nonlinear Differential Equations and their Applications 13[END_REF][START_REF] Fournais | Spectral Methods in Surface Superconductivity[END_REF][START_REF] Sandier | Progress in Nonlinear Differential Equations and their Applications[END_REF][START_REF] Sigal | Magnetic Vortices, Abrikosov Lattices and Automorphic Functions[END_REF]. Physically, this is because of the occurrence of mixed states where normal and superconducting regions may coexist. If the sample is a very long cylinder and the applied magnetic field is parallel to it, one may adopt a 2D description on a cross-section of the cylinder. One then distinguishes mainly two types of mixed phases:

• The vortex lattice (Abrikosov lattice [Abr]) where the normal regions take the form of vortices, and are arranged on a triangular lattice embedded in a sea of superconducting material; • The surface superconductivity state where the whole bulk of the sample is in the normal state, and superconductivity only survives close to the boundary. The second situation shall concern us here. In this regime, the magnetic field is very large, varies between two critical values H c2 and H c3 , and mostly penetrates the sample. That superconducting electrons may still exist because of boundary effects is a highly non-trivial fact, first derived in [SJdG]. At leading order, the phenomenon may be completely understood by considering the case of an infinite half-plane sample, with a straight boundary. It thus has some universal features: although superconducting electrons concentrate along the boundary, the geometry of the latter does not affect their distribution much. This is because the density of superconducting electrons essentially only varies in the direction normal to the boundary, as we proved rigorously in [START_REF] Correggi | On the Ginzburg-Landau Functional in the Surface Superconductivity Regime[END_REF][START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF], following several earlier contributions [START_REF] Almog | Nonlinear Surface Superconductivity in the Large κ Limit[END_REF][START_REF] Almog | The Distribution of Surface Superconductivity along the Boundary: on a Conjecture of X.B. Pan[END_REF][START_REF] Fournais | Superconductivity between Hc 2 and Hc 3[END_REF][START_REF] Fournais | Energy asymptotics for type II superconductors[END_REF][START_REF] Fournais | On the third critical field in Ginzburg-Landau theory[END_REF][START_REF] Fournais | Nucleation of bulk superconductivity close to critical magnetic field[END_REF][START_REF] Lu | Estimates of the upper critical field for the Ginzburg-Landau equations of superconductivity[END_REF][START_REF] Pan | Surface Superconductivity in Applied Magnetic Fields above Hc2[END_REF] (see [START_REF] Fournais | Spectral Methods in Surface Superconductivity[END_REF] for a review).

However, in order to prove some of the most refined results in [START_REF] Correggi | On the Ginzburg-Landau Functional in the Surface Superconductivity Regime[END_REF][START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF], we had to precisely estimate subleading order contributions to the energy. These do not share the universal character of the leading order, in that they depend on the sample, via the curvature of its boundary. This is reminiscent of earlier works, e.g., [START_REF] Fournais | Energy asymptotics for type II superconductors[END_REF] (see [START_REF] Fournais | Spectral Methods in Surface Superconductivity[END_REF]Chapters 13 and 15] for extensive references), where it has been proved that, when decreasing the magnetic field just below H c3 , superconductivity appears first where the curvature of the boundary is maximum. In this paper we aim at evaluating the effect of sample curvature in the whole regime of magnetic fields comprised between H c2 and H c3 . We shall give a simple expression of the curvature dependent contribution to surface superconductivity, which, again, appears only at subleading order in the energy, and thus requires rather refined estimates.

Our setting is the following: we consider the Ginzburg-Landau functional (in convenient units whose relation to other conventions is discussed 1 in [START_REF] Correggi | On the Ginzburg-Landau Functional in the Surface Superconductivity Regime[END_REF][START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF])

G GL ε [Ψ, A] = Ω dr ∇ + i A ε 2 Ψ 2 - 1 2bε 2 2|Ψ| 2 -|Ψ| 4 + 1 ε 4 |curlA -1| 2 . (1.1)
The domain Ω ⊂ R 2 represents the cross-section of an infinitely long cylinder of superconducting material. We assume that it is bounded, simply connected and that its boundary is smooth. The applied magnetic field is perpendicular to Ω.

We shall denote

E GL ε := min (Ψ,A)∈D GL G GL ε [Ψ, A], (1.2) with D GL := (Ψ, A) ∈ H 1 (Ω; C) × H 1 (Ω; R 2 ) , (1.3)
and denote by (Ψ GL , A GL ) a minimizing pair. We recall that |Ψ GL | 2 gives the local relative density of superconducting electrons (bound in Cooper pairs) and that curl A GL is the induced magnetic field in the sample. We are interested in the behavior of |Ψ GL | in the surface superconductivity regime

1 < b < Θ -1 0 , (1.4)
where Θ 0 is the minimal ground state energy of the shifted harmonic oscillator on the half-line:

Θ 0 := min α∈R min u 2 =1 +∞ 0 dt |∂ t u| 2 + (t + α) 2 |u| 2 . (1.5)
This corresponds to asking that the applied magnetic field varies between H c2 and H c3 , and we shall also assume that ε is a small parameter, in order to prove asymptotic results in the limit ε → 0. Physically this means we consider an "extreme" type-II superconductor.

In the parameter regime of our interest, the GL order parameter is concentrated near the boundary of the sample and the induced magnetic field is very close to the (constant) applied one. To leading order, superconducting electron pairs are uniformly distributed along the boundary as a function of the tangential variable. The main variations are in the direction normal to the boundary. Our main result in this note is an asymptotic estimate for |Ψ GL | 4 which exhibits a subleading (in ε) curvature-dependent correction. This shows that variations of the boundary's curvature influence the superconductivity distribution.

In our previous papers [START_REF] Correggi | On the Ginzburg-Landau Functional in the Surface Superconductivity Regime[END_REF][START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF], we have emphasized the role played by the simplified functional

E 1D k,α [f ] := c 0 | log ε| 0 dt(1 -εkt) |∂ t f | 2 + V k,α (t)f 2 -1 2b 2f 2 -f 4 (1.6)
where

V k,α (t) := (t + α -1 2 εkt 2 ) 2 (1 -εkt) 2 (1.7)
and c 0 is a (somewhat arbitrarily) fixed, large enough, constant. Note that in the limit ε → 0 the above expression reduces to a 1D nonlinear energy independent of the curvature, which is known to provide the leading order contribution to the GL asymptotics (see [START_REF] Correggi | On the Ginzburg-Landau Functional in the Surface Superconductivity Regime[END_REF] and references therein). This functional should be thought of as giving the GL energy to subleading order in the case of a sample Ω which is a "disc" of curvature k, i.e., either a disc of radius R = k -1 when k > 0 or the exterior of such a disc when k < 0. The variable t corresponds to the coordinate normal to the boundary (in units of ε). Denoting by s the tangential coordinate, f (t)e -iαs gives an ansatz for the GL order parameter in boundary coordinates. To get an optimal energy, this functional should be minimized with respect to both the function f and the number α, leading to an optimal profile f k , an optimal phase α(k) and an optimal energy E 1D ⋆ (k). Since we work in the regime ε → 0, it is natural to try to consider a perturbative expansion of E 1D ⋆ (k(s)). We then get that the leading order is given by the k = 0 functional (corresponding to a half-plane sample and extensively studied in the literature [AH, FHP, Pan]):

E 1D 0,α [f ] := +∞ 0 dt |∂ t f | 2 + (t + α) 2 f 2 -1 2b 2f 2 -f 4 , (1.8)
and the first correction by -εk times

E corr α [f ] := c 0 | log ε| 0 dt t |∂ t f | 2 + f 2 -α(t + α) - 1 b + 1 2b f 2 (1.9)
which is obtained by retaining only linear terms in εk when expanding (1.6). We shall denote E 1D 0 the minimum of E 1D 0,α [f ], and α 0 , f 0 a minimizing pair (α 0 is unique, f 0 also is, up to a sign). Note the mild abuse of notation: the k = 0 functional is well-defined even for c 0 = +∞, and we take this convention. Due to known decay estimates for f 0 (see, e.g., [START_REF] Correggi | On the Ginzburg-Landau Functional in the Surface Superconductivity Regime[END_REF]Proposition 3.3]), this creates only an exponentially small discrepancy in ε, provided c 0 is a large enough constant.

We previously proved energy estimates relating the full GL energy to the infimum of the above 1D, curvature-dependent functional. Since our method was local, and the GL energy density is related to |Ψ GL | 4 , it is natural to expect a result about the distribution of the latter quantity from the energy estimate. The goal of this paper is to provide this estimate.

Let us first introduce scaled boundary coordinates: the surface superconductivity layer Ãε := {r ∈ Ω | τ ≤ c 0 ε| log ε|} , (1.10) where τ := dist(r, ∂Ω), (1.11) can be mapped to

A ε := {(s, t) ∈ [0, |∂Ω|] × [0, c 0 | log ε|]} (1.12) via a diffeomorphism Φ.
For technical reasons (see Remark 1.2 below), we can only evaluate D |Ψ GL | 4 with the desired precision in the case that the set D looks "rectangular" in boundary coordinates. Let then D ⊂ Ω be a measurable set independent of ε such that

Φ(D ∩ Ãε ) = [s D , s ′ D ] × [0, c 0 | log ε|] (1.13) for some s D , s ′ D ∈ [0, |∂Ω|].
Notice that this implies that the boundary of D intersects ∂Ω with π/2 angles. Our main result is the following: Theorem 1.1 (Curvature dependence of the order parameter). Let Ψ GL be a GL minimizer and D ⊂ Ω be a measurable set such that (1.13) holds. Denote s → k(s) the curvature of ∂Ω as a (smooth) function of the tangential coordinate s. For any 1 < b < Θ -1 0 , in the limit ε → 0,

D dr |Ψ GL | 4 = ε C 1 (b)|∂Ω ∩ ∂D| + ε 2 C 2 (b) ∂D∩∂Ω ds k(s) + o(ε 2 ), (1.14)
where ds stands for the 1D Lebesgue measure along ∂Ω and

C 1 (b) = -2bE 1D 0 = +∞ 0 dt f 4 0 > 0 (1.15) C 2 (b) = 2b E corr α 0 [f 0 ] = 2 3 bf 2 0 (0) -2bα 0 E 1D 0 .
(1.16)

Moreover, for |b -Θ -1 0 | small enough (independently of ε), C 2 (b) > 0. The leading order term in (1.14) had previously been computed [FK, Kac, Pan], with a less explicit expression of the constant C 1 (b) however.

Remark 1.1. (Concentration of Cooper pairs.) For obvious physical reasons it would be preferable to have an estimate of D |Ψ GL | 2 , which would directly give information on the distribution of Cooper pairs close to the boundary of the sample. Unfortunately, our method, which is mostly energy-based, does not provide this. An estimate on D |Ψ GL | 4 is easier to obtain because more directly linked to the concentration of the energy density

e GL ε (r) := ∇ + i A GL ε 2 Ψ GL 2 - 1 2bε 2 2|Ψ GL | 2 -|Ψ GL | 4 + b ε 4 curlA GL -1 2 , (1.17)
as we shall see below. The above theorem still indicates that, to leading order, |Ψ GL | is concentrated evenly along the boundary of Ω, and that the first correction to this effect is directly proportional to the curvature function k(s). Notice that, in order to exploit (1.17), the use of the variational equation solved by the GL minimizer is crucial and therefore our result does not extend directly to low energy configurations, although one would expect it.

Remark 1.2. (Convergence as measures.)

It would be natural to reformulate the result in terms of convergence of measures, i.e., by stating that 1 ε

1 ε |Ψ GL | 4 dr -C 1 (b)ds(r) -→ ε→0 C 2 (b)k(s)ds(r), (1.18)
in the sense of measure, where ds(r) is the 1D Lebesgue measure along the boundary of Ω. However, due to the restriction (1.13) on the shape of the set D, the statement (1.14) is weaker than (1.18). We in fact expect (1.14) to be wrong as stated if the set D intersects the boundary with an angle = π/2. Indeed, the result is based on integrating the energy density on lines normal to the boundary that cover the full extent of the physical region. These integrals give the energy of the 1D model (1.6) that enters the main formula. One can easily decompose D ∩ A ε into the union of sets of the form (1.13) and some remainder which is more triangular-like. In the latter regions we cannot follow the proof procedure to obtain a simple expression. This does not affect the leading order of the result, as noted in [FK, Pan], since the energy contained in the triangular regions is small relatively to the leading order. It is however not small compared to the correction we isolate in (1.14) because the area of the triangular regions can easily be seen to be O(ε 2 ). Integrating |Ψ GL | 4 on such an area gives a O(ε 2 ) contribution and thus ruins the result.

Technically speaking, Assumption (1.13) enters the proof in the estimate of the boundary integral appearing in Lemma 4.1. It allows to exploit a new pointwise estimate of the tangential derivative of |Ψ GL | 2 that we prove in Lemma 4.2. It is important to remark that such an estimate does not hold for the normal component of the gradient (see the discussion preceding Lemma 4.2).

Remark 1.3. (Limiting regimes.) In [START_REF] Fournais | Nucleation of bulk superconductivity close to critical magnetic field[END_REF]Corollary 1.3] a similar result about the boundary behavior of surface superconductivity is proven for magnetic fields slightly below the second critical one H c2 , i.e., for b → 1 -, b ≤ 1. In this estimate the leading order is the same as in (1.14) but the first correction is proportional to the area of the set |D| and is due to the bulk behavior of the superconductor.

The regime where b → Θ -1 0 at the same time as ε → 0 was studied in details in [START_REF] Fournais | Energy asymptotics for type II superconductors[END_REF] (see in particular [START_REF] Fournais | Energy asymptotics for type II superconductors[END_REF]Theorem 1.4] and the discussion thereafter). The behavior of the GL functional becomes approximately linear in this limit, and therefore the whole machinery of spectral theory of linear operators can be exploited to extract a lot of details about the boundary behavior. In particular, according to the asymptotics of Θ -1 0 -b when ε → 0, superconductivity can be either uniformly distributed all over the boundary or concentrated close to the points of maximal curvature. The first order correction to the boundary behavior is however not known, although thanks to the approximate linearity of the problem, estimates can be formulated in terms of the integral of

|Ψ GL | 2 instead of |Ψ GL | 4 .
Remark 1.4. (Sign of the curvature correction) Unfortunately we are not able to determine the sign of the correction in (1.14). Based on the results of [START_REF] Fournais | Energy asymptotics for type II superconductors[END_REF] we have just recalled, we conjecture that

E corr α 0 [f 0 ] > 0 for any 1 < b < Θ -1
0 . This would mean that points with large curvature attract more superconductivity. We can prove this conjecture only when b is close (independently of ε) to Θ -1 0 , see Lemma 2.3.

Note that the 2D setting we consider here corresponds to an infinite 3D cylinder with a magnetic field parallel to the axis. In a more general 3D setting, the angle between the magnetic field and the surface of the sample also plays a role in the distribution of surface superconductivity, see [FKP] and references therein.

The rest of the paper contains the proof of Theorem 1.1. In Section 2 we first discuss the perturbative expansion of the 1D ground state energy, obtain the expression (1.16) and prove that it is positive close to the third critical field. Section 3 contains a result of the same form as Theorem 1.1 where |Ψ GL | 4 is replaced by the GL energy density (1.17). Finally, in Section 4, which is the more involved of the three, we deduce our main result from the estimate of the energy density.

Acknowledgments. M.C. acknowledges the support of MIUR through the FIR grant 2013 "Condensed Matter in Mathematical Physics (Cond-Math)" (code RBFR13WAET). N.R. acknowledges the support of the ANR project Mathostaq (ANR-13-JS01-0005-01).

We are indebted to one of the anonymous referees of our previous paper [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF], whose remarks motivated the present investigation.

Correction to the 1D Energy

The asymptotic expansion of the 1D energy that is behind (1.14) is a direct consequence of first-order perturbation theory. We first prove that the functional (1.9) gives the subleading order of the 1D energy.

Let us recall the notation: for any k ∈ R, E 1D ⋆ (k) is the infimum of (1.6) w.r.t. both f and α, while f 0 and α 0 stand for the minimizing pair of (1.8).

Lemma 2.1 (Perturbative expansion of the 1D energy).

As ε → 0

E 1D ⋆ (k) = E 1D 0 -εkE corr α 0 [f 0 ] + O(ε 3/2 | log ε| γ ) = - 1 2b +∞ 0 f 4 0 -εkE corr α 0 [f 0 ] + O(ε 3/2 | log ε| γ ) (2.1)
for some fixed γ > 0 where E corr is defined as in (1.9).

Proof. We first take f 0 , α 0 as a trial pair for the functional (1.6), which yields

E 1D ⋆ (k) ≤ E 1D 0 -εkE corr α 0 [f 0 ] + O(ε 2 )
by expanding E 1D k,α 0 [f 0 ] and using that α 0 , f 0 do not depend on ε. In order to estimate the remainders it suffices to use the exponential decay of

f 0 [CR1, Proposition 3.3], which implies that ∞ 0 dt t n f 2 0 (t) = O(1),
for any n ∈ N. The decay estimate also implies that we make no significant error by considering the k = 0 functional as being defined on the whole half-line, provided c 0 is large enough.

Next we write

E 1D ⋆ (k) = E 1D k,α(k) [f k ] = E 1D 0,α(k) [f k ] -εkE corr α(k) [f k ] + O(ε 2 ) ≥ E 1D 0 -εkE corr α(k) [f k ] + O(ε 2 ),
where we use the variational principle defining E 1D 0 . Finally it easily follows from the estimates of [CR2, Proposition 1] that

E corr α(k) [f k ] = E corr α 0 [f 0 ] + O(ε 1/2 | log ε| γ )
for some γ > 0. Gathering the previous inequalities, we get an upper and a lower bound to E 1D ⋆ (k) which together give (2.1). That

E 1D 0 = - 1 2b +∞ 0 dt f 4 0 (2.2)
follows by mutliplying the variational equation

-∂ 2 t f 0 + (t + α 0 ) 2 f 0 = 1 b (1 -f 2 0 )f 0 (2.3)
by f 0 and integrating.

The expression (1.16) is somewhat unwieldy, but can be simplified a lot:

Lemma 2.2 (Simple formula for the 1D energy correction).

Let α 0 , f 0 be a minimizing pair for the 1D functional (1.8) at k = 0. Then, for all

1 < b < Θ -1 0 , E corr α 0 [f 0 ] = 2 3 (1 -α 2 0 b) + α 0 2b +∞ 0 dt f 4 0 = 1 3 f 2 0 (0) -α 0 E 1D 0 .
(2.4)

Proof. Let us first recall two useful identities:

+∞ 0 dt (t + α 0 )f 2 0 (t) = 0 (2.5)
expresses the optimality of α 0 , see, e.g., [START_REF] Fournais | Superconductivity between Hc 2 and Hc 3[END_REF]Eq. (3.20)] or [CR1, Lemma 3.1], while

+∞ 0 dt t(t + α 0 )f 2 0 (t) = +∞ 0 dt |∂ t f 0 | 2 + 1 4b f 4 0 (2.6)
is a virial identity, equivalent to [START_REF] Fournais | Superconductivity between Hc 2 and Hc 3[END_REF]Eq. (3.22)]. It is obtained by noting that

E α 0 (ℓ) := inf f E 1D 0 1 √ ℓ f (•/ℓ) = inf f E 1D 0 [f ]
for all ℓ and thus

∂ ℓ E α 0 (ℓ)| ℓ=1 = 0.
Using the Feynman-Hellmann principle to evaluate the latter expression one gets (2.6).

We now start the computation:

Step 1. We claim that

E corr α 0 [f 0 ] = +∞ 0 dt t b - 3α 0 4b f 4 0 (t) + +∞ 0 dt 2t 3 - 2t b -2α 2 0 t + α 0 b f 2 0 (t). (2.7)
First, using the virial identity (2.6) we obtain

E corr α 0 [f 0 ] = +∞ 0 dt (t -α 0 )|∂ t f 0 | 2 + t 2b - α 0 4b f 4 0 (t) - t b f 2 0 .
Next, integrating by parts and using the Neumann boundary condition for

f 0 , +∞ 0 dt (t -α 0 )|∂ t f 0 | 2 = - +∞ 0 dt (t -α 0 ) 2 ∂ t f 0 ∂ 2 t f 0 .
Inserting the variational equation ( 2.3) and integrating by parts again we deduce

+∞ 0 dt (t -α 0 )|∂ t f 0 | 2 = +∞ 0 dt (t -α 0 ) 2 (t + α 0 ) + (t + α 0 ) 2 (t -α 0 ) - t -α 0 b + t -α 0 2b f 2 0 (t) f 2 0 (t) + α 4 0 2 f 2 0 (0) - α 2 0 2b f 2 0 (0) + α 2 0 4b f 4 0 (0).
The claim then follows from the identity (see, e.g., [CR1, Proof of Lemma 3.3])

f 2 0 (0) = 2 -2α 2 0 b (2.8)
and a bit of algebra.

Step 2. Next we compute that

+∞ 0 dt (t + α 0 ) 3 f 2 0 = - 1 2b +∞ 0 dt (t + α 0 )f 4 0 + 1 3 1 -α 2 0 b , (2.9) following [FH3, Proof of Lemma 3.2.7]. Let us denote H α 0 := -∂ 2 t + (t + α 0 ) 2 + 1 b f 2 0 -1 b
and recall that f 0 is a ground state for this Schrödinger operator, with eigenvalue 0. For any function v we therefore have, integrating by parts

f 0 , H α 0 v = v ′ (0)f 0 (0).

If we apply this with

v = (t + α 0 ) 2 f ′ 0 -(t + α 0 )f 0 , we find f 0 , H α 0 v = α 4 0 + α 2 0 b (f 2 0 (0) -1) -1 f 2 0 (0).
On the other hand, using the variational equation,

-∂ 2 t v = -3(t + α 0 )f ′′ 0 -(t + α 0 ) 2 f ′′′ 0 = -(t + α 0 ) 5(t + α 0 ) 2 + 3 b (f 2 0 -1) f 0 -(t + α 0 ) 2 (t + α 0 ) 2 + 3 b f 2 0 -1 b f ′ 0 , so that H α 0 v = -6(t + α 0 ) 3 f 0 -4 b (t + α 0 )(f 2 0 -1)f 0 -2 b (t + α 0 ) 2 f 2 0 f ′ 0 .
Multiplying this by f 0 and integrating yields

f 0 , H α 0 v = -6 +∞ 0 dt (t + α 0 ) 3 f 2 0 - 3 b +∞ 0 dt (t + α 0 )f 4 0 + α 2 0 2b f 4 0 (0)
where we use also (2.5) and an integration by parts. Equating the two different expressions of f 0 , H α 0 v we obtained and noting that, since

f 2 0 (0) = 2 -2α 2 0 b, f 2 0 (0) -α 4 0 - α 2 0 b (f 2 0 (0) -1) + 1 + α 2 0 2b f 4 0 (0) = f 2 0 (0) = 2 1 -α 2 0 b
yields (2.9).

Step 3. From (2.9), (2.5) and (2.6) we get

+∞ 0 dt t 3 f 2 0 = +∞ 0 dt (t + α 0 ) 3 -3α 2 0 t -3α 0 t 2 -α 3 0 f 2 0 = 1 6 f 2 0 (0) + +∞ 0 dt - 1 2b (t + α 0 )f 4 0 -3α 0 |∂ t f 0 | 2 - 3α 0 4b f 4 0 -α 3 0 f 2 0 .
Inserting in (2.7) and using (2.5) again we find

E corr α 0 [f 0 ] = 1 3 f 2 0 (0) + +∞ 0 dt - 13α 0 4b f 4 0 -6α 0 |∂ t f 0 | 2 + α 0 b - 2t b -2α 2 0 t -2α 3 0 f 2 0 = 1 3 f 2 0 (0) + +∞ 0 dt - 13α 0 4b f 4 0 -6α 0 |∂ t f 0 | 2 + 3 α 0 b f 2 0 .
(2.10)

Finally we note that

+∞ 0 dt |∂ t f 0 | 2 = E 1D 0 - 1 2b +∞ 0 dt f 4 0 + +∞ 0 dt 1 b -(t + α 0 ) 2 f 2 0 = - 1 b +∞ 0 dt f 4 0 + +∞ 0 dt 1 b -(t + α 0 ) 2 f 2 0 whereas (2.5) and (2.6) together imply +∞ 0 dt (t + α 0 ) 2 f 2 0 = +∞ 0 dt |∂ t f 0 | 2 + 1 4b f 4 0 so that, combining the two identities, +∞ 0 dt |∂ t f 0 | 2 = +∞ 0 dt 1 2b f 2 0 - 5 8b f 4 0 .
Inserting this in (2.10) and recalling once more (2.8) and (2.2) this yields the final expressions (2.4).

Unfortunately we are not able to determine the sign of the energy correction from the expressions we found. However, when b → Θ -1 0 we have more information: it is known that f 2 0 scales as (1 -bΘ 0 ) 1/2 . It then immediately follows from Lemma 2.2 that the correction must be positive for b close enough to Θ -1 0 : Lemma 2.3 (Sign of the correction close to H c3 ). There exists

1 < b 0 < Θ -1 0 such that, for all b 0 < b < Θ -1 0 , E corr α 0 [f 0 ] > 0.
(2.11)

Proof. We recall from [FH3, Section 3.2] that the minimum in (1.5) is achieved by a unique pair u 0 , α opt with

α opt = - √ Θ 0 and u 0 normalized in L 2 . At b > Θ -1 0 it is easy to see that f 0 ≡ 0. When b → Θ -1
0 one should therefore expect f 0 → 0, in which case the quartic term becomes a second order correction. One should thus expect that the solution is close to that of the linear problem, the only effect of the quartic term being to fix the overall normalization. More precisely, following the techniques of [FH3, Section 14.2.2] one can show that α 0 -→ b→Θ -1 0 -Θ 0 and that

u 0 4 4 1 -bΘ 0 1/2 f 0 -→ b→Θ -1 0 u 0 .
It is easy to see that the latter convergence holds in the quadratic form domain of the harmonic oscillator. Using standard elliptic estimates, one can upgrade this to any Sobolev or Hölder norm. In particular, convergence holds in L ∞ and L 4 , so that

f 2 0 (0) = (1 -bΘ 0 ) u 0 4 4 u 2 0 (0)(1 + o b (1))
where we denoted by o b (1) a quantity going to 0 as b → Θ -1 0 , and

+∞ 0 dt f 4 0 = (1 -bΘ 0 ) 2 u 0 4 4 (1 + o b (1)).
From Lemma 2.2 we have that

E corr α 0 [f 0 ] = 1 3 f 2 0 (0) + α 0 2b +∞ 0 dt f 4 0 = 1 -bΘ 0 3 u 0 4 4 u 2 0 (0) + 3α 0 2b (1 -bΘ 0 ) = 1 -bΘ 0 3 u 0 4 4 u 2 0 (0) + o b (1) > 0 since u 0 (0) > 0 is independent of b.

Estimates of the Energy Density

Our main estimate on the order parameter is obtained by exploiting the link between |Ψ GL | 4 and the GL energy density. We discuss first the asymptotics for the latter.

Proposition 3.1 (Estimates for the energy density). Let e GL ε be the GL energy density defined in (1.17) and D ⊂ Ω be a measurable set satisfying (1.13). Under the same assumptions and with the same notation as in Theorem 1.1 we have, as ε → 0

D dr e GL ε = ε -1 E 1D 0 |∂Ω ∩ ∂D| -E corr α 0 [f 0 ] ∂D∩∂Ω ds k(s) + O(ε 1/2 | log ε| γ ). (3.1)
Proof. This is an adaptation of the method developed in [START_REF] Correggi | On the Ginzburg-Landau Functional in the Surface Superconductivity Regime[END_REF][START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF]. First we recall from [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF]Lemma 4] the energy lower bound

E GL ε ≥ 1 ε G Aε [ψ] -Cε 2 | log ε| 2 , (3.2)
where ψ is, up to a phase factor, the GL order parameter in boundary coordinates and the reduced functional is

G Aε [ψ] := |∂Ω| 0 ds c 0 | log ε| 0 dt (1 -εk(s)t) |∂ t ψ| 2 + 1 (1 -εk(s)t) 2 |(ε∂ s + ia ε (s, t)) ψ| 2 - 1 2b 2|ψ| 2 -|ψ| 4 (3.3) with a ε (s, t) := -t + 1 2 εk(s)t 2 + εδ ε , (3.4)
and

δ ε := γ 0 ε 2 - γ 0 ε 2 , γ 0 := 1 |∂Ω| Ω dr curl A GL , (3.5) ⌊ • ⌋ standing for the integer part.
This lower bound may in fact (with identical proof) be localized, yielding

ε D dr e GL ε ≥ s ′ D s D ds c 0 | log ε| 0 dt (1 -εk(s)t) |∂ t ψ| 2 + 1 (1 -εk(s)t) 2 |(ε∂ s + ia ε (s, t)) ψ| 2 - 1 2b 2|ψ| 2 -|ψ| 4 -Cε 3 | log ε| 2 . (3.6) We next split the interval [s D , s ′ D ] into N ε = O(ε -1 ) sub-intervals [s n , s n+1 ], n = 1, . . . , N ε , of side length O(ε). The convention is that s D = s 1 and s ′ D = s Nε+1 . This gives a decomposition of [s D , s ′ D ] × [0, c 0 | log ε|] into N ε rectangular cells C n , n = 1, . . . , N ε . Arguing as in [CR2, Lemma 6] we then deduce ε D dr e GL ε ≥ s ′ D s D ds E 1D ⋆ (k(s)) + Nε n=1 E n [u n ] -Cε 2 | log ε| ∞ (3.7)
where, within the n-th cell,

ψ(s, t) =: u n (s, t)f n (t) exp -i αn ε + δ ε s , (3.8)
and the reduced functionals E n are defined as

E n [u] := Cn dsdt (1 -εk n t) f 2 n |∂ t u| 2 + 1 (1-εknt) 2 |ε∂ s u| 2 -2εb n (t)J s [u] + 1 2b f 2 n 1 -|u| 2 2 , (3.9) with b n (t) := t + α n -1 2 εk n t 2 (1 -εk n t) 2 , (3.10) and J s [u] := (iu, ∂ s u) = i 2 (u * ∂ s u -u∂ s u * ) . (3.11)
The mean curvature in the n-th cell is denoted k n , with f n and α n the minimizing profile and phase for the associated functional (1.6). Inserting the result of Lemma 2.1 into (3.7) we find

D dr e GL ε ≥ ε -1 E 1D 0 |∂Ω ∩ ∂D| -E corr α 0 [f 0 ] ∂D∩∂Ω ds k(s) + ε -1 Nε n=1 E n [u n ] -Cε| log ε| ∞ . (3.12)
We next adapt the strategy of [CR2, Proof of Lemma 7] to estimate the reduced functionals from below:

Nε n=1 E n [u n ] ≥ -Cε 3/2 | log ε| γ , (3.13)
thereby concluding the proof of the lower bound corresponding to (3.1). This is a long procedure that we will not recall in details. We shall emphasize the only point that has to be modified, due to the fact that we now bound from below the energy density in the set D instead of the full GL energy.

Step 1 of the proof of [CR2, Lemma 7] adapts with no modification, leading to

E n [u n ] ≥ ε tn,ε 0 dtF n (t) [J t [u n ](s n+1 , t) -J t [u n ](s n , t)] + d ε Cn dsdt (1 -εk n t) f 2 n |∂ t u n | 2 + 1 (1-εknt) 2 |ε∂ s u n | 2 + O(ε ∞ ), (3.14)
for some d ε ∼ | log ε| -4 , denoting

J t [u] := (iu, ∂ t u) = i 2 (u * ∂ t u -u∂ t u * ) and F n (t) = 2 t 0 dη (1 -εk n η)b n (η)f 2 n (η).
We then follow Step 2 of the same proof to combine and estimate the boundary terms produced by the use of Stokes' formula in Step 1 (terms on the first line of the above formula). In this procedure it is crucial to sum boundary terms living on the same cell boundary. Since here s 1 = s Nε+1 there is obviously a need for a different estimate of the terms located on the corresponding boundaries, i.e., those of the original set D. This is the only point where we depart slightly from the method of [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF] and rely on more refined inequalities.

We proceed as follows (say for the n = 1 term, located on the boundary corresponding to s 1 = s D ): let χ be a smooth cut-off function depending only on s with

χ(s 1 ) = 1, |χ| ≤ 1, supp(χ) ⊂ C 1 , |∇χ| ≤ Cε -1 .
Since our cells have side-length O(ε) in the s direction, the last two requirements are obviously compatible. Intergrating by parts in the s variable we get

t1,ε 0 dt F 1 (t)J t [u 1 ](s 1 , t) = t1,ε 0 dt χ(s 1 )F 1 (t)J t [u 1 ](s 1 , t) = C 1 dsdt χF 1 ∂ s J t [u 1 ] + C 1 dsdt F 1 ∂ s χJ t [u 1 ].
(3.15)

We drop the subscripts 1 for shortness. To handle the first term in the above we note that

∂ s J t [u] = i 2 (∂ s u∂ t u * -∂ s u * ∂ t u) + i 2 (u∂ s ∂ t u * -u * ∂ s ∂ t u)
and hence a further integration by parts in t yields

C dsdt χF ∂ s J t [u] = - i 2 C dsdt χ∂ t F (u∂ s u * -u * ∂ s u).
Note that the boundary terms vanish by definition of F . We also have

|∂ t F | ≤ C| log ε|f 2 and thus ε C dsdt χF ∂ s J t [u] ≤ Cε| log ε| C dsdt f 2 |u||∂ s u| ≤ Cδ| log ε| C dsdt f 2 |u| 2 + Cδ -1 | log ε| C dsdt f 2 |ε∂ s u| 2 ≤ Cδε| log ε| 2 + Cδ -1 ε 2 | log ε| γ ≤ Cε 3/2 | log ε| γ
where we use that f 2 |u| 2 = |ψ| 2 ≤ 1 plus the fact that |C| = O(ε| log ε|), recall the estimate [CR2, Eq. (6.15)] and have chosen δ = ε 1/2 | log ε| γ for the final optimization.

For the second term in (3.15) we write, using essentially the same ingredients (in particular [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF]Eq. (6.15)

]), ε C dsdt F ∂ s χJ t [u] ≤ C C dsdt f 2 |u||∂ t u| ≤ Cδ C dsdt f 2 |u| 2 + Cδ -1 C dsdt f 2 |∂ t u| 2 ≤ Cδε| log ε| + Cδ -1 ε 2 | log ε| γ ≤ Cε 3/2 | log ε| γ .
Combining the previous estimates, we obtain

ε t1,ε 0 dt F (t)J t [u](s 1 , t) = O(ε 3/2 | log ε| γ )
and a similar estimate for the term located on the boundary s = s Nε+1 . Dealing with the other boundary terms as in [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF] concludes the proof of (3.13).

At this stage we have the lower bound corresponding to (3.1),

D dr e GL ε ≥ ε -1 E 1D 0 |∂Ω ∩ ∂D| -E corr α 0 [f 0 ] ∂D∩∂Ω ds k(s) + O(ε| log ε| γ ) (3.16)
and by the same method also

D c dr e GL ε ≥ ε -1 E 1D 0 |∂Ω ∩ ∂D c | -E corr α 0 [f 0 ] ∂D c ∩∂Ω ds k(s) + O(ε 1/2 | log ε| γ ). (3.17)
On the other hand, combining the energy estimate of [CR2, Theorem 1] and Lemma 2.1 we have the global estimate

Ω dr e GL ε = ε -1 E 1D 0 |∂Ω| -E corr α 0 [f 0 ] ∂Ω ds k(s) + O(ε| log ε| γ ).
Combining with (3.17) we deduce

D dr e GL ε = Ω dr e GL ε - D c dr e GL ε ≤ ε -1 E 1D 0 |∂Ω ∩ ∂D| -E corr α 0 [f 0 ] ∂D∩∂Ω ds k(s) + O(ε 1/2 | log ε| γ )
which we combine with (3.16) to complete the proof.

From the Energy Density to the Order Parameter

We now conclude the proof of Theorem 1.1 by adding the following to Proposition 3.1:

Proposition 4.1 (Energy density versus order parameter).

Under the assumptions and with the notation of Theorem 1.1 and Proposition 3.1, we have

D dr e GL ε = - 1 2bε 2 D dr |Ψ GL | 4 + o(1). (4.1)
The proof is split in two lemmas. First we have a general result which does not require the set under consideration to be rectangular: Let S ⊂ Ω be a measurable set. Then, with the previous notation

S dr e GL ε + 1 2bε 2 S dr |Ψ GL | 4 = 1 2 ∂S dσ ∇|Ψ GL | 2 • ν + o(1) (4.2)
with ν the outward-pointing normal to ∂S.

Proof. We first note that

1 ε 4 S dr curlA GL -1 2 = O(ε| log ε| 3 ). (4.3)
Indeed, using the elliptic estimate (see [START_REF] Fournais | Spectral Methods in Surface Superconductivity[END_REF]Eq. (11.51)])

curlA GL -1 C 1 (Ω) = O(ε)
and the fact that curlA GL = 1 on ∂Ω, we deduce that in the full boundary layer (1.10) we have curlA

GL -1 = O(ε 2 | log ε|)
and thus

S∩ Ãε dr curlA GL -1 2 = O(ε 5 | log ε| 3 ).
The part of the integral located in S ∩ Ãc ε is of much lower order, as follows from the usual Agmon estimates, for instance [START_REF] Fournais | Spectral Methods in Surface Superconductivity[END_REF]Eq. (12.10)], and we deduce (4.3).

At the level of precision we aim at we may thus neglect the magnetic kinetic energy:

S dr e GL ε = S dr ∇ + i A GL ε 2 Ψ GL 2 - 1 2bε 2 2|Ψ GL | 2 -|Ψ GL | 4
+ O(ε| log ε| 3 ).

Next we recall that since Ψ GL is a minimizer for G GL ε we have the first Ginzburg-Landau variational equation:

-∇ + i A GL ε 2 2 Ψ GL + 1 bε 2 Ψ GL |Ψ GL | 2 -1 = 0. (4.4)
Combining with the identity

1 2 ∆|Ψ GL | 2 = ℜ(Ψ GL ∆Ψ GL ) + |∇Ψ GL | 2 we deduce 1 2 ∆|Ψ GL | 2 = ∇ + i A GL ε 2 Ψ GL 2 + 1 bε 2 |Ψ GL | 2 |Ψ GL | 2 -1 . (4.5)
Integrating over S we obtain S dr e GL ε = -

1 2bε 2 S |Ψ GL | 4 + 1 2 S dr ∆|Ψ GL | 2 (4.6)
and the proof is complete since of course

S dr ∆|Ψ GL | 2 = ∂S dσ ∇|Ψ GL | 2 • ν.
Applying the previous lemma with D = S, our task should now be to estimate the boundary term in the right-hand side of (4.2). It is similar to terms showing up in [FK, Proof of Lemma 6.1], but using the estimates therein shows at best that it is of order O(1), a remainder that we cannot afford. A technical novelty in the present paper is thus to show that this term is in fact o(1), provided the set D is rectangular in boundary coordinates.

We certainly have

∂D dσ ∇|Ψ GL | 2 • ν = (∂D)∩Aε dσ ∇|Ψ GL | 2 • ν + O(ε ∞ )
by the usual decay estimates, where we recall that A ε is defined in (1.12). Splitting the curve L D := Φ(∂D ∩ A ε ) (which is a rectangle in boundary coordinates) in a part L s D parallel to the boundary of ∂Ω and a part L t D normal to the boundary of ∂Ω we have

∂D dσ ∇|Ψ GL | 2 • ν = L t D dσ ∇|Ψ GL | 2 • ν + O(ε ∞ ). ( 4
∇Ψ GL + i A GL ε 2 • ν = 0, on ∂Ω
satisfied by Ψ GL . The other part of L s D is deep in the region where the order parameter decays exponentially and may thus be neglected. The new key ingredient is that on L t D we can prove

∇|Ψ GL | 2 • ν = ∂ s |Ψ GL | 2 ≤ Cε a-1 (4.8)
for some a > 0. This is natural in view of the results of [START_REF] Correggi | On the Ginzburg-Landau Functional in the Surface Superconductivity Regime[END_REF][START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF]: the variations in the s-direction should be much smaller than those in the t-direction, which happen on a scale ε. Combining this with the fact that the length of this part of the boundary is O(ε| log ε|) we get the desired estimate. We are in fact not able to prove (4.8) in all the boundary region A ε and will thus have to split the line integral on L t D into two pieces. Let us introduce the following subset of A ε where the estimate (4.8) is proven in the next Lemma 4.2:

Ã> := r ∈ Ω f 0 (τ /ε) ≥ ε 1/6 .
(4.9)

Any power of ε strictly smaller than 1/4 would do the job but we fix it equal to 1/6 for concreteness. Recall that τ = dist(r, ∂Ω). As before we denote by A > = Φ( Ã> ), the set Ã> in boundary coordinates and it is easy to see that

A > = [0, ∂Ω] × [0, t > ] , t > ≫ 1. (4.10)
In fact, exploiting the available pointwise bounds on f 0 (see for example [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF]Eq. (A.6)]), one immediately verifies that

t > ≥ 1 √ 3 | log ε|(1 + o(1)). (4.11)
Before stating the pointwise estimate mentioned above, let us stress that a similar bound cannot hold for the normal component of the gradient of Ψ GL : the estimate

∂ t |Ψ GL | 2 ∝ ε -1
is optimal. The different behavior of the s and t derivatives will be apparent in the proof of the following Lemma.

Lemma 4.2 (Estimate of the tangential derivative).

As ε → 0 ∂ s Ψ GL (Φ(s, τ )) 2 L ∞ (A>) = O(ε -5/6 | log ε| ∞ ).
(4.12)

Proof. The starting point is the variational equation (4.4) satisfied by Ψ GL : setting as in [CR2, Sect. 5] ψ(s, t) = Ψ GL (Φ(s, εt))e -iφε(s,t) , (4.13) where the explicit expression of the gauge phase is given in [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF]Eq. (5.4)], one gets (4.14) i.e., thanks to the choice of the gauge the magnetic field is now purely tangential. The explicit expression of à can be easily recovered in terms of φ ε (see, e.g., [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF]Eq. (5.6)]), but the most important point is the estimate Ã(s, t)

-∂ 2 t ψ + 1 (1-εk(s)t) 2 -iε∂ s + Ã ε 2 ψ = 1 b 1 -|ψ| 2 ψ,
+ εt L ∞ (Aε) = O(ε 2 | log ε| 2 ),
which follows from a priori estimates on A GL as [START_REF] Correggi | On the Ginzburg-Landau Functional in the Surface Superconductivity Regime[END_REF]Eq. (4.23)]. In addition the explicit expression [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF]Eq. (5.6)] also implies that ∂ s à ∞ = O(ε| log ε|). Notice that here we are exploiting the smoothness of ∂Ω and the fact that the curvature is infinitely differentiable. Plugging the ansatz ψ(s, t) = f 0 (t)e -i α 0 ε s u(s, t), (4.15)

for some unknown function u and with f 0 and α 0 the minimizing density and phase of the half-plane functional (1.8), we get

-∂ 2 t (f 0 u) + f 0
(1-εk(s)t) 2 (-iε∂ s + α 0 + t + O(ε| log ε|)) 2 u = 1 b 1 -f 2 0 |u| 2 f 0 u. Exploiting now the variational equation (2.3) for f 0 and dividing by f 0 > 0, we obtain where we have used the upper bound |u| ≤ f -1 0 ≤ ε -1/6 and the estimate |1 -|u|| = O(ε 1/4 | log ε| ∞ ) in A > which follows from [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF]Theorem 2]. Rescaling now also the s variable by setting ξ = s/ε and denoting v(ξ, t) = u(εξ, t), we can apply the Gagliardo-Nirenberg inequality [START_REF] Nirenberg | On Elliptic Partial Differential Equations[END_REF]p. 125]

-∂ 2 t u -2 f ′ 0 f 0 ∂ t u - 1 (1-εk(s)t) ε 2 ∂ 2 s u -2iε(α 0 + t)∂ s u =
∇ ξ,t v ∞ ≤ C ∆ ξ,t v 1/2 ∞ 1 -|v| 1/2 ∞ + 1 -|v| ∞ ,
which implies ∇v ∞ = O(ε 1/6 | log ε| ∞ ) and therefore

(∂ t + ε∂ s ) |u| L ∞ (A>) ≤ ∇v ∞ = O(ε 1/6 | log ε| ∞ ).
(4.16)

The result on Ψ GL then follows from the identities (4.13) and (4.15).

Putting together the results of Lemma 4.1 and Lemma 4.2, we are now in position to complete the proof of the main result of this Section:

Proof of Proposition 4.1. The estimate (4.2) and the properties of the set D yield (4.7) and therefore 1), (4.17) since

S dr e GL ε + 1 2bε 2 S dr |Ψ GL | 4 = 1 2 ∂S dσ ∇|Ψ GL | 2 • ν + o(1) = L t D dσ ∇|Ψ GL | 2 • ν + o(1) = L t D ∩ Ã> dσ ∇|Ψ GL | 2 • ν + o(
L t D ∩ Ãc > dσ ∇|Ψ GL | 2 • ν = 2 Φ(L t D ∩ Ãc > ) dτ |ψ(s, τ /ε)|∂ s |ψ(s, τ /ε)| = 2 Φ(L t D ∩ Ãc > ) dτ |ψ(s, τ /ε)|f 0 (τ /ε)∂ s |u(s, τ /ε)| ≤ Cε -1 ψ L ∞ (A c > ) Φ(L t D ∩ Ãc > )
where we have used [START_REF] Correggi | Boundary behavior of the Ginzburg-Landau order parameter in the surface superconductivity regime[END_REF]Eq. (6.2)]. Now the Agmon estimate for ψ stated, e.g., in [CR2, Eq. (5.5)] yields

ψ L ∞ (A c > ) ≤ Ce -At> ≤ C exp -1 √ 3
| log ε| ≪ | log ε| -1 , thanks to (4.11). Hence we conclude that

L t D ∩ Ãc > dσ ∇|Ψ GL | 2 • ν ≤ C| log ε| ψ L ∞ (A c > ) = o(1),
and (4.17) is proven. For the rest of the boundary integral it suffices to apply (4.12):

L t D ∩ Ã> dσ ∇|Ψ GL | 2 • ν = Φ(L t D )∩A> dσ ∂ s Ψ GL (Φ(s, τ )) 2 = O(ε 1/6 | log ε| ∞ ).

  Lemma 4.1 (Reduction to a boundary term).

  .7) Indeed, on the part of L s D which coincides with ∂Ω we have ∇|Ψ GL | • ν = 0 by taking the real part of the Neumann boundary condition

  |u| 2 + O(ε| log ε| 2 ) u. Since (see [CR1, Lemma A.1]) f ′ 0 f 0 = O(| log ε| 3 ),inside A > the above equation yields the estimate∂ 2 t + ε 2 ∂ 2 s u ≤ C | log ε| 3 |(∂ t + ε∂ s ) u| + |1 -|u|| |u| ≤ C | log ε| 3 |(∂ t , ε∂ s )u| + ε 1/12 | log ε| b


Hélène Perrier

Jérémy Levallois

David Coeurjolly

Jean-Philippe Farrugia

Jean-Claude Iehl

Jacques-Olivier Lachaud

Interactive Curvature Tensor Visualization on Digital Surfaces

Keywords: Isosurface Visualization, Digital Geometry, Curvature Estimation, GPU

Interactive visualization is a very convenient tool to explore complex scientific data or to try different parameter settings for a given processing algorithm. In this article, we present a tool to efficiently analyze the curvature tensor on the boundary of potentially large and dynamic digital objects (mean and Gaussian curvatures, principal curvatures, principal directions and normal vector field). More precisely, we combine a fully parallel pipeline on GPU to extract an adaptive triangulated isosurface of the digital object, with a curvature tensor estimation at each surface point based on integral invariants. Integral invariants being parametrized by a given ball radius, our proposal allows to explore interactively different radii and thus select the appropriate scale at which the computation is performed and visualized.

Introduction

Volumetric objects are being more and more popular in many applications ranging from object modeling and rendering in Computer Graphics to geometry processing in Medical Imaging or Material Sciences. When considering large volumetric data, interactive visualization of those objects (or isosurfaces) is a complex problem. Such issues become even more difficult when dynamic volumetric datasets are considered. Beside visualization, we are also interested in performing geometry processing on the digital object and to explore different parameter settings of the geometry processing tool. Here, we focus on curvature tensor estimation (mean/Gaussian curvature, principal curvatures directions. . . ). Most curvature estimators require a parameter fixing the scale at which the computation is performed. For short, such parameter (integration radius, convolution kernel size. . . ) specifies a scale for analyzing the object surface, and is naturally related to the amount allowed perturbations on input data. As a consequence, when using such estimators, exploring different values of this parameter is mandatory. However, this is usually an offline process, due to the amount of computations that need to be done.

Contributions

In this work, we propose a framework to perform interactive visualization of complex 3D digital structures combined with a dynamic curvature tensor estimation. We define a fully data parallel process on the GPU (Graphics Processor Unit) to both efficiently extract adaptive isosurface and compute per vertex curvature tensor using Integral Invariants estimators. This system allows us to visualize curvature tensor in real-time on large dynamic objects. Our approach combines a GPU implementation of pointerless octrees to represent the data, an adaptive viewpoint-dependent mesh extraction, and a GPU implementation of integral invariant curvature estimators.

Related works Extracting and visualizing isosurface on volumetric data has been widely investigated since the seminal Marching Cubes approach by Lorensen and Cline [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF]. This approach being data parallel, GPU implementation of this method is very efficient [START_REF] Tatarchuk | Real-time isosurface extraction using the gpu programmable geometry pipeline[END_REF]. However, such technique generates a lot of triangles which is not well suited to large digital data. Hence, adaptive approaches have been proposed in order to optimize the triangulated mesh resolution according to the object geometry or the viewpoint. In this topic, many solutions have been developed in Computer Graphics [START_REF] Shu | Adaptive marching cubes[END_REF][START_REF] Schaefer | Dual marching cubes: Primal contouring of dual grids[END_REF][START_REF] Lengyel | Voxel-based terrain for real-time virtual simulations[END_REF][START_REF] Lewiner | Fast generation of pointerless octree duals[END_REF][START_REF] Lobello | Out-of-core adaptive iso-surface extraction from binary volume data[END_REF]. The method developed by Lengyel et al. [START_REF] Lengyel | Voxel-based terrain for real-time virtual simulations[END_REF] suits best our needs. It combines an octree space partitioning with new Marching Cubes configurations to generate adaptive meshes on CPU from static or near static data. We propose here a full GPU pipeline inspired by this method, that maintains a view dependent triangulation. Our high framerate allows us to inspect dynamic 3D data in real-time.

Curvature estimation on discrete or digital surface has also been widely investigated. In [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF], authors propose digital versions of Integral Invariant estimators [START_REF] Pottmann | Principal curvatures from the integral invariant viewpoint[END_REF][START_REF] Pottmann | Integral invariants for robust geometry processing[END_REF] in order to estimate the complete curvature tensor (mean/Gaussian curvatures, principal curvatures, principal curvature directions, normal vector field) on digital surfaces. Such approaches are based on an integration principle using a ball kernel of a given radius. Additionally, authors have demonstrated that these estimators have multigrid convergence properties. In this article, we also propose an efficient GPU implementation of such estimators to visualize such curvature fields in real-time and to interactively change the value of the kernel radius.

In this paper, we first present the previous works on curvature tensor estimation. Then, we propose an efficient GPU approach to extract a triangulated isosurface from a digital object. Finally, we propose a fully-parallel GPU pipeline to compute curvature tensor in real-time and present our results.

× x BR (x) ∂X X h • z ∈ X ∩ BR (x) BR (x) × x h ∂X ∂[Z] h

Curvature Tensor Estimation

In our context, we consider digital shapes (any subset Z of Z d ) and boundaries of digital shapes Bd(Z). We denote by G h (X) the Gauss digitization of a shape

X ⊂ R d in a d-dimensional grid with grid step h, i.e. G h (X) := {z ∈ Z d , h • z ∈ X}. For such digitized set Z := G h (X), [Z] h is a subset of R d corresponding
to the union of hypercubes centered at h • z for z ∈ Z with edge length h. By doing so, both ∂X and ∂[Z] h are topological boundaries of objects lying in the same space (see Fig. 1-b). Note that the combinatorial digital boundary Bd(Z) of Z made of cells in a cellular Cartesian complex (pointels, linels, surfels, . . . ), can be trivially embedded into R d such that it coincides with ∂[Z] h .

In [START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF], authors define a 2D digital curvature estimator κR and a 3D digital mean curvature estimator ĤR based on the digital volume estimator Vol(Y, h) := h d Card(Y ) (area estimator Area(Y, h) in dimension 2): Definition 1 Given the Gauss digitization Z := G h (X) of a shape X ⊂ R 2 (or R 3 for the 3D mean curvature estimator), digital curvature estimators are defined for any point x ∈ R 2 (or R 3 ) as:

∀0 < h < R, κR (Z, x, h) := 3π 2R - 3 Area(B R/h (x/h) ∩ Z, h) R 3 , (1)
ĤR (Z, x, h) := 8 3R - 4 Vol(B R/h (x/h) ∩ Z, h) πR 4 . ( 2
)
where B R/h (x/h) is the ball of digital radius R/h centered on digital point (x/h) ∈ Z.

Such curvature estimators have multigrid convergence properties [START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF]: when the digital object becomes finer and finer, i.e. when the digitization step h tends to zero, the estimated quantities on ∂[G h (X)] h converges (theoretically and experimentally) to the associated one on ∂X in O h

1 3
for convex shapes with at least C 3 -boundary and bounded curvature (setting R := kh 1 3 for some k ∈ R). In [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF], authors have also defined 3D digital principal curvature estimators κR Definition 2 Given the Gauss digitization Z := G h (X) of a shape X ⊂ R 3 , 3D digital principal curvature estimators are defined for any point x ∈ R 3 as:

∀0 < h < R, κR 1 (Z, x, h) := 6 πR 6 ( λ2 -3 λ1 ) + 8 5R , (3)
κR 2 (Z, x, h) := 6 πR 6 ( λ1 -3 λ2 ) + 8 5R , (4)
where λ1 and λ2 are the two greatest eigenvalues of the covariance matrix of B R/h (x/h) ∩ Z.

The covariance matrix needs to compute digital moments of order 0, 1 and 2 (see Equation 20 of [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF] for more details). These estimators are proven convergent in O h

1 3
when setting the ball radius h in R = kh 1 3 , where k is a constant related to the maximal curvature of the shape, on convex shapes with at least C 3 -boundary and bounded curvature [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF]. Additionally, eigenvectors associated to λ1 and λ2 of the covariance matrix are principal curvature direction estimators ŵR 1 and ŵR 2 . The smallest eigenvector corresponds to the normal direction nR at x. Convergence results can be found in [START_REF] Lachaud | Robust and convergent curvature and normal estimators with digital integral invariants[END_REF].

It has been shown that the radius of the ball depends on the geometry of the underlying shape. In [START_REF] Levallois | Parameter-free and multigrid convergent digital curvature estimators[END_REF], a proposal was made for a parameter-free estimation of the radius of the ball by analyzing the shape w.r.t. the local shape geometry using maximal digital straight segments of the digital boundary. In [START_REF] Levallois | Scale-space feature extraction on digital surfaces[END_REF], these estimators have been analyzed in scale-space (for a range of radii) for a given digital shape. This allows to detect features of the shape thanks to the behavior of estimators on singularities. As a consequence, for all these integral invariant based approaches, we need to consider different ball radius which could be time consuming when implemented on CPU. We propose here a fully parallel implementation on GPU allowing us to change the radius and thus update the estimated quantities in real-time.

Isosurface Extraction on GPU

In this section, we detail the adaptive isosurface extraction algorithm. The proposed approach uses an octree representation of the input object on which an adaptive Marching Cube builds the isosurface efficiently. Such hierarchical representation of the object allows us to handle large datasets and to locally adapt the level of details w.r.t. the geometry or camera position. We first present the octree representation and then the isosurface extraction.

Linear Octree Representation

Representing a hierarchical structure on GPU is usually challenging since such a data parallel component is unable to handle recursivity. Efficient spatial tree encoding can be achieved using pointerless structures such as linear quadtrees or octrees Gargantini [4]. This structure indexes each cell by a Morton code: the code of children cells are defined by the code of the parent suffixed by two bits (in dimension 2, three bits in dimension 3) (see Figure 2-left). A cell's code encodes its position w.r.t. its parent cell and its complete path to the tree root. Hence, the tree is fully represented as a linear vector of its leaves. Furthermore, cell operations such as subdivision, merging can be efficiently implemented using bitwise operations on the Morton code. In the following, we use the GPU friendly implementation proposed by Dupuy et al. [START_REF] Dupuy | GPU Pro 5, chap. Quadtrees on the GPU[END_REF].

Data parallel and adaptive mesh generation

Using this spatial data structure, a triangulated mesh can be constructed using Marching Cubes [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF] (MC for short): the triangulation is generated from local triangle patches computed on local cell configurations. Such approach is fully parallel and easy to implement on the GPU. However, since adjacent cells may not have the same depth in the octree, original Lorensen and Cline's rules need to be updated (see Figure 2-right). Many authors have addressed this problem both for primal and dual meshes [START_REF] Shu | Adaptive marching cubes[END_REF][START_REF] Schaefer | Dual marching cubes: Primal contouring of dual grids[END_REF][START_REF] Lengyel | Voxel-based terrain for real-time virtual simulations[END_REF][START_REF] Lewiner | Fast generation of pointerless octree duals[END_REF][START_REF] Lobello | Out-of-core adaptive iso-surface extraction from binary volume data[END_REF].

In the following, we use the algorithm proposed by Lengyel et al. [START_REF] Lengyel | Voxel-based terrain for real-time virtual simulations[END_REF]. First, this approach forces the octree structure to make sure that the depth difference between any two adjacent cells is at most one. Then, Lengyel et al. introduce the concept of transition cells. Those cells are defined to be inserted between two neighboring octree cells of different depth. With specific MC configurations for their triangulation, a crack free mesh can be extracted.

Similarly to original MC algorithm, this approach is well suited to a GPU implementation: given a set of cells (a vector of morton codes), each triangle patch can be extracted in parallel for both regular cells and transition cells.

Level of Details Criteria and Temporal Updates

As illustrated in Figure 2-middle, we propose a viewpoint dependent criterion to decide if a cell needs to be refined: the closer we are to the camera, the finer the cells are. Such a criterion is also well suited to GPU since it can be evaluated independently on every cell. Figure 3 illustrates our level of details (LoD for short) criterion. In dimension 2, if α denotes the viewing angle, an object at a distance d from the camera has a projected size on screen of 2 • d • tan(α). (see Figure 3-left). Our distance criterion is based on the ratio (visibility ratio in the following) between the cell diameter l(c) (power of 2 depending on the depth), and its projected size. For a given cell c, split and merge decision are based on this visibility ratio:

c is split if its children cells have a visibility ratio greater than constant k ; c and its sibling cells are merged if their parent cell c has a visibility ratio lower than k ; otherwise, the cell c stays for the next frame. Using such a criterion, split and merge decisions are computed in parallel from the morton codes of all cells.

Once decisions have been made, a new set of cells is sent to the mesh generation step described above. Finally, before constructing the triangulation from remaining cells and transition cells, geometrical culling is performed in order to skip the triangle patch construction for cells that are not visible. Figure 4 illustrates the overall fully data parallel pipeline. We first present design principles for the GPU implementation and then present our approaches. The general idea is to perform an Integral Invariant computation on GPU at each vertex of the generated triangulated mesh. Since GPU have massively parallel architectures, we can do all those computations in parallel to obtain a very efficient implementation. Please note that when the LoD criterion is removed, each MC vertex is exactly centered at a surfel center. At different depth of the octree, MC vertices still correspond to surfel centers of subsampled versions of the input object. Hence, Integral Invariant framework defined in Section 2 is consistent with the triangulated surface obtained on GPU: triangles are used for visualization purposes but all computations are performed on the digital object G h (X) for estimators (1), ( 2), ( 3) and (4).

To implement the integration on B R (x) ∩ X (Fig. 1), several strategies have been evaluated.

Per Vertex Real-time Computation on GPU

Naive Approach. A first simple solution consists in scanning all digital points, at a given resolution, lying inside the integration domain (Fig. 5-lef t) and then estimating the geometrical moment as the sum of the geometrical moments of each elementary cubes lying in the intersection (see Eq. ( 1) to ( 4)). This is exactly similar to what is done on the CPU. On the GPU, we can exploit mipmap textures to obtain multi-resolution information. If the input binary object is stored in a 3D texture, GPU hardware constructs a multi-resolution pyramid (mipmap) such that 8 neighboring voxel intensities at level l are averaged to define the voxel value at level l + 1. If the level 0 corresponds to the binary input object, at a given level l, a texture probe at a point (x, y, z) returns the fraction of G h (X) belonging to the cube of center (x, y, z) and edge length 2 l . As a consequence, we can approximate the volume of B R (x) ∩ X by considering mipmap values at a given resolution l (Fig. 5-right). In this case, errors only occurs for cells lying inside X (with density 1) not entirely covered by B R (x). Furthermore, the mipmap texture can be used to design, using a single texture probe, a fast inclusion test of a given cell c at level l into the shapes: we say that c is in X if its density (retrieved from the texture probe at level l) is greater than 1/2. The idea here is to mimic a kind of adaptive Gauss digitization process.

Hierarchical Decomposition. Using the hierarchical nature of the mipmap texture, we could also consider hierarchical decompositions of B R (x). The idea is to decompose the ball into mipmap cells of different resolution in order to limit the number of texture access (important bottleneck on GPU hardwares) and to get better approximated quantities. On the GPU, computing a hierarchical decomposition of B R (x) at each point x is highly inefficient since the hardware optimizes the parallelism only when the micro-program described in the shader has a predictive execution flow. Hence, implementing the recursive algorithm (or its de-recursified version) requires a lot of branches (conditional Fig. 5. Here are the three approaches we used to evaluate the integrand. (left) is the naive approach, where we sum all the digital points that are inside the integrand domain (middle) shows the hierarchical decomposition of the ball, thus limiting the number of texture probes to do. (right) illustrates the dynamic approach, we probe the textures at a higher resolution and we refine it at each frame until we reach the finer level.

if tests) in the flow. We have thus considered a fast multi-resolution approximation as illustrated in Fig. 5-middle: for a given radius R, we precompute a hierarchical decomposition of B R . Such decomposition is made of octree cells at different resolutions. At a given MC vertex x, the algorithm becomes simple since we just scan the ball cells and estimate the area (or other moments) from the mipmap values associated to this cell. Note that these precomputed cells in the B R decomposition may not be aligned with mipmap cells. However, we can use GPU which interpolates mipmap values if we probe at non-discrete positions. For a given radius R, the expected number of cells in B R is in O(log R). Implementation details on the hierarchical octree representation can be found in the supplementary material (see Section 6).

Dynamic Refinement. In this last approach, we want to optimize the interactivity and the execution flow or parallelism on GPU. The integration is simply computed using a regular grid at different resolution l (Fig. 5-right). The shader code becomes trivial (simple triple loop) and a lot of texture fetches are involved, but interactivity can be easily obtained. Indeed, we consider the following multipass approach:

1. When the surface geometry has been computed, we set the current level l at an high level l := l max of the mipmap pyramid. 2. We compute the integrals (and the curvature tensor) at the mipmap level l and send the estimated quantities to the visualization stage. 3. If the user changes the camera position or the computation parameters, we return to step 1.

4. If the current framerate is above a given threshold, we decrease l and return to step 2 if the final l 0 level has not been reached yet.

A consequence of the multi-pass approach is that when there is no interaction (camera settings, parameters), the GPU automatically refines the estimation.

Even if step 2 is quite expensive, in O R 2 l
3 , the interactive control in step 3 considerably improves the interactive exploration of the tensor with fast preliminary approximations which are quickly refined.

Compared to the hierarchical approach, no precomputation is required for a given radius R. As a consequence, we could even locally adapt the ball radius to the geometry (for instance following the octree depth of the current MC vertex). In the next section, we evaluate the performances of both approaches. We first evaluate curvature estimations on a full-resolution geometry obtained by disabling the LoD criterion at the mesh generation step. Figure 6-top shows results of curvature tensor estimation (mean, Gaussian, first and second principal directions) on "OctaFlower" at level l 0 (considered as our ground truth in the following). Figure 6-bottom shows mean curvature estimation in real-time on a dynamic object. For this case, we simply evaluate the implicit expression at each vertex of a cell to decide if such cell is included or not into B R (x) ∩ X instead of updating the mipmap texture of densities at each time step. For all illustrations, we use directly normal vectors computed by algorithm discussed in Section 2.

Experiments

Full resolution experiment

Then, we compare the approximations made by computing curvature from level l ≥ l 0 in Figure 7. We can see that results using approximation seems to quickly converge to ground truth results. Table 1 shows numerical results of L ∞ error (maximal absolute difference for all vertices) and L 2 error (mean squared errors of all vertices), as well the number of mipmap texture fetches. We can also see that the number of texel fetch, required to compute the curvature, reduces drastically when computing approximations. However, at higher levels, approximation errors become more important. Those levels are thus never used in practice, they are presented here to illustrate how the refining process converges to the l 0 level.

We also compare them with the hierarchical algorithm (as discussed in Section 4.1). This method introduces a higher error when compared with l 0 . This is due to the precomputation of the subdivision that no longer ensures to fetch data at the center of a mipmap cell. The GPU interpolation reduces the bias, but it does not remove it.

Fully adaptive evaluation

When dealing with large datasets, we cannot expect real-time curvature tensor estimation if we consider the full resolution geometry, due to the huge amount of data to process. By computing a dynamically refined approximate curvature tensor estimation joined with an adaptive triangulation (as discussed in Section 4.1), we manage to maintain a real-time framerate by giving control over the amount of data to process at each frame. In Figure 8, we compare timings (in logarithmic scale) for a triangulation that is dynamically refined according to its distance to the camera. We measured those timings using the multiresolution regular and the hierarchical algorithm.

First we can note that the required time to compute the ground truth curvature for an object is usually as high as the time required to extract its geometry and to compute all of the above levels. Using approximations is thus mandatory. Another advantage with approximations is that it allows us to get a visualization of our object as soon as we run the application. This allows for real-time interactions with the object, required in order to change the visualized quantity, curvature radius, etc.

It is also visible in Figure 8 that a hierarchical decomposition greatly reduces the curvature computation time, especially with big radii. However, due to the precomputation and the current hierarchical structure, this algorithm is biased and creates an error (presented in Table 1) that needs to be considered.

Figure 9 shows curvature computation and exploration in real-time on large . Timings in milliseconds (in logscale) obtained while visualizing an adaptive triangulation on three objects -"OctaFlower" (digital domain of 130 3 ), "Snow microstructures" (233 3 ) and "XYZ-Dragon" (510 3 ) -by computing the curvature with a regular grid (left) and with hierarchical algorithm (right), with two different radii for each object. In orange color : time required to extract the triangulation. In blue color : time required to compute the curvature tensor at different levels: from l4 (light blue) to l0 (dark blue). Timings are given using a NVIDIA GeForce GTX 850M GPU.

datasets: "XYZ-Dragon" (with digital domain of 512 3 ) and "Snow microstructures" (233 3 ).

Conclusion and Discussion

In this article, we have proposed a fully data parallel framework on GPU hardware which combines an adaptive isosurface construction from digital data with a curvature tensor estimation at each vertex. Using this approach, we can explore in real-time different curvature measurements (mean, Gaussian, principal directions) with different ball radii on potentially large dynamic dataset. Our proposal relies on both a linear octree representation with Morton codes and an efficient integral computation on GPU. The source code and additional ma-terial (video, . . . ) are available on the project website (https://github.com/ dcoeurjo/ICTV).

Fig. 1 .

 1 Fig. 1. Integral invariant computation (left) and notations (right) in dimension 2 [2].

1 and κR 2

 12 on Z ⊂ Z 3 based on digital moments:

Fig. 2 .

 2 Fig. 2. Morton codes associated to cells of a linear quatree. Each morton code of a child cell is obtained by adding a suffix to its parent code (left). The adaptive representation consists of quadtree cells whose depth is view point dependent (middle). Finally, adaptive Marching Cubes is used to generate the triangulation (right).

Fig. 3 .

 3 Fig. 3. Notations (left) and adaptive meshing in dimension 2 using the LoD distance and angular criterion (right).

Fig. 4 .

 4 Fig. 4. GPU pipeline summary. Data buffers are represented in green and computations in red. Each computation retrieves data from a buffer and fills a new one.

Fig. 6 .

 6 Fig. 6. First row : Mean and Gaussian curvature estimation, (zoom of) first and second principal directions estimation on "OctaFlower" with a digital domain of 130 3 . Second row: Mean curvature computed in real-time (around 20 FPS) on a dynamic object.

Fig. 7 .

 7 Fig. 7. Illustration of mean curvature computation on "OctaFlower" (digital domain of 130 3 ) using mipmap approximation with different levels: l3, l2, l1 and l0 (i.e. no approximation).

Fig. 8

 8 Fig.8. Timings in milliseconds (in logscale) obtained while visualizing an adaptive triangulation on three objects -"OctaFlower" (digital domain of 130 3 ), "Snow microstructures" (233 3 ) and "XYZ-Dragon" (510 3 ) -by computing the curvature with a regular grid (left) and with hierarchical algorithm (right), with two different radii for each object. In orange color : time required to extract the triangulation. In blue color : time required to compute the curvature tensor at different levels: from l4 (light blue) to l0 (dark blue). Timings are given using a NVIDIA GeForce GTX 850M GPU.

Fig. 9 .

 9 Fig. 9. Left column: Mean curvature, first and second principal directions and normal vector field estimation on "Snow microstructures" (233 3 and R = 8). Right column: Mean curvature, first and second principal directions and normal vector field estimation on "XYZ-Dragon" (510 3 and R = 8). Normal vectors are colored with a mapping of their component to RGB color space.

Table 1 .

 1 Comparison of number of texture fetches, L2 and L∞ error obtained on "OctaFlower" with a digital domain of 130 3 when computing mean curvature with two radii: 8 and 16, with hierarchical algorithm (H) and l ≥ l0 approximation algorithms. The object is triangulated at full resolution with a regular grid and contains 282,396 vertices.

			H	l4	l3	l2	l1	l0
	Number of	R = 8	1468		2	28	260 2104
	texture fetches	R = 16	5706	2	28	90	2120 17080
	L∞ error	R = 8	0.051		0.306	0.085	0.047	0
	(w.r.t. l0)	R = 16	0.053	0.146	0.041	0.017	0.005	0
	L2 error	R = 8 3.51e-05		2.67e-04 7.57e-05 2.02e-05	0
	(w.r.t. l0)	R = 16 4.43e-05 1.39e-04 4.16e-05 1.57e-05 2.07e-06	0

3

) using mipmap approximation with different levels: l3, l2, l1 and l0 (i.e. no approximation).

This work has been mainly funded by DigitalSnow ANR-11-BS02-009, KIDICO ANR-2010-BLAN-0205 and PRIMES Labex ANR-11-LABX-0063/ ANR-11-IDEX-0007 research grants.


Abdelbaki Djouambi
email: djouambi_abdelbaki@yahoo.fr

Abdelfatah Charef
email: afcharef@yahoo.com

Alina Voda Besançon

APPROXIMATION AND SYNTHESIS OF NON INTEGER ORDER SYSTEMS

Keywords: Fractional systems, Approximation, Identification and Distribution of Relaxation Times

published or not. The documents may come

Approximation and synthesis of non integer order systems

Abdelbaki Djouambi, Abdelfatah Charef, Alina Voda

INTRODUCTION

Fractional order systems or fractal systems are defined in terms of the shape of their spectra which is proportional to reciprocal frequency and which is generally called 1/f β spectra [START_REF] Mandelbrot | Some noises with 1!f spectrum, a bridge between direct current and white noise[END_REF] and [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF]. Such spectra have been observed in a tremendous variety of dissimilar natural phenomena.

Numerous studies have been made about this phenomenon ever since it was pointed out by (Van Der Ziel, 1950). Many empirical and mathematical models have also been introduced to describe these processes with 1/f β type spectra. During the past decade, we have seen tremendous growth of study in this subject especially toward the applications in meteorology, astronomy, geology, physiology and many other fields. Recently, fractal analytical techniques and methodologies have been used by scientists in different fields in order to explain the ubiquitous 1/f β processes which have been characterizing a wide range of natural phenomena, see [START_REF] Dutta | Low frequency fluctuations in solids: 1/f noise[END_REF] and [START_REF] Sun | A unified approach to represent metal electrode interface[END_REF]. These scientists have found that these processes are mediated by the fractal nature of the phenomena themselves. Hence, these observations suggest a link between fractals and system theory.

The equation that is commonly used to describe some phenomenon has been expressed by an inverse power law equation as:

m s K s F = ) ( (1)
Where, k is a positive constant, s = jw, and m is a positive real number.

Hence many of the complex fractal systems which have been represented as the inverse power law equation can best be represented by the Cole-Cole function and described in the Laplace transform domain by the explicit transfer function [START_REF] Cole | Dispersion and Absorption in Dielectrics, Alternating Current Characterization[END_REF]:

a s k s F m + = ) ( exp (2)
Where, k and a is a positive constant, jw s =

, and m is a positive real number. This type of expression gives a much suitable mathematical representation to the natural phenomena so that its low frequency magnitude is finite instead of infinite. This is usually the case for the magnitude of natural phenomena at very low frequency.

The analysis of the response to elementary input of these fractional order systems shows that they posses a remarkable performance as compared to the integer order systems, see [START_REF] Oustaloup | Frequencyband complex non integer differentiator: Characterization and synthesis[END_REF]) [START_REF] Oustaloup | From fractal robustness to non integer approach in vibration insulation[END_REF], [START_REF] Charef | Fractal systems as represented by singularity function[END_REF], [START_REF] Sun | Fractal system -a time domain approach[END_REF] and [START_REF] Charef | Fractional Order Feedback Control Systems, JIEEEC2001 4 th conference[END_REF].

In this paper we present a new method to approximate the transfer function given by equation (2) using the distribution of relaxation times function.

The materials collected in this paper are grouped into two main parts: Part I summarizes the approximation by a rational function of equation ( 2) which we found to be the basic element for representation of fractional order systems. Part II deals with the synthesis with the Foster's network of the approximated fractional system.

APPROXIMATION

Many of real physical systems can be batter represented by an explicit fractional transfer function given by the following transfer function [START_REF] Cole | Dispersion and Absorption in Dielectrics, Alternating Current Characterization[END_REF]:

m s s F ) ( 1 1 ) ( 0 τ + = (3)
where

ℜ ∈ m .
Because of their mathematical representation in the frequency domain are irrational function, direct analysis methods and corresponding time domain behavior seems difficult to handle. For the purpose of identification, analysis, synthesis and simulation of such systems, the need arises for a rational function approximation. Two cases were presented here:

Case 1: 0 < m <1
To extend Debye model [START_REF] Debye | Polar Molecules[END_REF]) which has only one relaxation time, the distribution of relaxation times has also been derived. This concept was introduced by Fuross and Kirkwood [START_REF] Fuross | Electrical Properties of Solids VIII-Dipole Moments in polyvinyl Chloride Biphenyl Systems[END_REF] derived a formula to obtain the distribution of relaxation times function ) (τ G directly from the original transfer function as fellow:

τ τ τ d s G s F ∫ ∞ + = 0 1 ) ( ) ( (4)
( ) [ ] ( ) [ ] π τ τ π π τ m m m G - - ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = 1 cos log cosh 1 sin 2 1 ) ( 0 , (5)
By sampling the distribution function G(τ) given in equation ( 4) for a limited bandwidth at logarithmically equidistant points i τ , for i=1, …, 2N-1 we can write:

( ) ( ) ∑ - = - = 1 2 1 ) ( N i i i G G τ τ δ τ τ (6)
where, ( ) ( ) and 1 > λ is the ratio of time constant to the next one or the ratio of a pole to the previous one, where the poles are:

i N i - = λ τ τ 0 , for i=1, 2,…, 2N-1,
i i s τ / 1 = , for i=1, 2,… , 2N-1.
Then, the transfer function of equation ( 3) can be approximated by:

( ) ∑ + = ≈ - = 1 2 1 1 ) ( ) ( N i i i N s G s F s F τ τ (7)
Using MATLAB Optimization Toolbox, The values of the ratio λ are calculated according to the following criterion: 200 4.190 4.190 As an illustration example, we have chosen a fractional order system which is represented by the following function: Using direct method time domain simulation, see [START_REF] Oustaloup | La dérivation non entière, théorie, synthèse et application[END_REF], the step response of the original system and that of the approximated one were presented in figure (3), the figure shows that the two curves are superimposed. As an illustrative example, we choose a fractional order system which is represented by the following function: ; the figures show that the two curves are overlapping with a small difference in phase responses due to the approximation error. 7) can be easily replaced by a respective driving point impedance of the following RC circuit type: The results obtained in this paper validate the link between non integer derivation and the distributed parameter systems. To approximate the Cole-Cole function by a first-order time-varying differential equation whose solution is also the inverse Laplace-Transform of the original fractional system, the distribution of relaxation times was used. By doing this transformation, the singularity function method for the explicit fractional system was used. The circuit representation consisted in a cascade of RC circuits of the singularity function.

[ ] { }⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - = ∈ ) ( ) (
i i i i i i C sR R s Z I V + = = 1 ) (
Z 1 (s) Z 2 (s) Z 2N-1 (s) I V R i C i Vi I i

  The figures (1) and (2) show the Bode plots of the original function and that of its rational approximation function with N=30 and 718 . 2 = λ . The figures show that the two curves are overlapping.

Fig

  Fig.1. Magnitude plot of ) ( ) ( s andF s F N

Fig. 3

 3 Fig.3Step response of the considered system and its approximation

Figure

  Figure (4) and figure (5) show the Bode plot of the original function and that of its approximation with 79 . 1 = e and 30 = N; the figures show that the two curves are overlapping with a small difference in phase responses due to the approximation error.

Fig

  Fig.4. Magnitude plot of ) ( ) ( s andF s F a

Fig

  Fig. 6 Elementary RC circuit

  of N and m . The values of this ratio for four different values of m and four different values of N are given in table (Tab.1).

	max min , w max w in the frequency band of interest min arg λ F N w	jw w	∈	jw 6 10 , 10 -F [	6	]	rd/s
	for different values Tab.1: The optimal values of the ratio λ , the maximum
	magnitude error (	eg max ) and the maximum phase error
	eph max for different values of m and N, with	τ	0 =	1	.
		N		10	20		30	40
	m	λ		10.434 3.040		2.650 2.711
	0.2	max	eg	8.851	2.776		0.200 0.026
		max	eph	3.880	4.197		0.890 0.100
		λ		7.126	2.716		2.718 2.719
	0.4	max	eg	6.520	0.395		0.010 0.003
		max	eph	4.260	4.031		0.040 0.006
		λ		4.727	2.717		2.721 2.721
	0.7	max	eg	4.630	0.023		0.010 0.010
		max	eph	9.170	0.340		0.080 0.090
		λ		4.250	3.844		3.844 3.844
	0.9	max	eg	2.370	1.550		1.550 1.550
		max	eph	5.650	4.

Proceedings of the 2nd IFAC Workshop on Fractional Differentiation and its Applications Porto,Portugal, July 19-21, 2006