M Siroux

E Tang-Kwor

D Especel

S Matteï

Emissivity measurements by using sine wave thermal modulation. Radiometric and calorimetric approaches

This paper presents a new technique for measuring total emissivities at room temperature involving a multifrequency sine wave thermal modulation. The infrared output signal as well as the sample's front side temperature are measured. Experimental data and a one dimensional model are compared, yielding identified values of diffusivity a, Biot number Bi and E, the latter being defined as the ratio of the sample to reference emissivity. Using Bi, emissivity can be also reached, so emissivity of a thermally thick sample can be reached through two different approaches. Results have been obtained for a 5mm thick sample of PVC, and it is shown that they are well concordant with integrated directional measurements.

Introduction

When modeling heat transfer the total radiative properties of surfaces have to be known. These properties are highly dependent on the surface conditions. They can not be found easily in the literature, and must be measured in working conditions. Emission measurement techniques may be separated into two groups : radiometric measurements and calorimetric measurements [START_REF] Sacadura | Proc. 9 th Int. Heat Transfer Conf[END_REF]. Steady state radiometric techniques are not convenient for low temperature emissivity measurements, because of the sample reflected flux. Thermal periodic techniques have been used mainly for measuring thermal diffusivity (Cowan, 1961), [START_REF] Ranchy | High Temp[END_REF], however a periodic radiometric method has been developed for emissivity measurements at room temperature (Matteï et al. 1994). This method is limited to thin sample with known surface temperature. The knowledge of surface temperature is needed for transient or calorimetric emission measurements too. This paper presents a periodic technique for measuring total emissivities of thick materials at low temperature, without emissive surface temperature measurement. This technique allows to reach total emissivity through two different approaches : radiometry and calorimetry.

Thermal behavior of a modulated strip shaped sample

Consider a strip shaped sample (Fig. 1) the front side temperature of which, T FS , is modulated :

T t T T f t FS ( ) cos( ) = + + 0 0 2π ϕ (1)
We define the thermal transfer function of the strip ~( ) H f 0 by the relation :

~( ) ~( ) ~( ) H f T f T f BS FS 0 0 0 = (2)
where ~( ) T f BS 0 and ~( ) T f FS 0 are the complex temperature of respectively the back side and the front side of the sample.

H is a complex number, the modulus of which characterizes the temperature attenuation, and the phase, the temperature phase shift between the two sides of the strip. When ~( ) H f is equal to one then we say that the sample is thermally thin. For a modulated temperature one dimensional conduction model, the heat balance equation can be written : H q e Bi sh q e q e Bi ch q e

∂ ∂ π 2 2 2 0 ~T x j f a T - = (
(5)

e is the sample thickness. Bi is the Biot number given by Bi h e k .

. h is the heat transfer coefficient on the sample back side which is the sum of a convective coefficient and a radiative coefficient. k is the conductivity. q is given by q j f a = 2π

, where a is the sample diffusivity.

Application to emissivity measurements

Principle of the method

Consider a total infrared detector facing the back side of the sample the temperature of which is in the form given by Eq. 1. The measured radiative flux is the sum of an emitted flux and a reflected one. At temperature modulation frequency f 0 , the reflected flux is not modulated. If ~( ) K f is the transfer function of the IR detection system (including geometrical parameters), the IR output signal f 0 component is given by :

~( ) ~( ). . ~( ) U f K f s T f BS 0 0 4 0 = ε σ π (6)
ε s is the sample back side emissivity. Its geometrical properties depend on the IR detection system optical path. We note ~( ) R f 0 the ratio :

~( ) ~( ) ~( ) R f U f T f FS 0 0 4 0 = (7)
If temperature amplitudes are small, ~( ) T f BS 4 0 can be linearised, and if front side and back side mean temperature are equal, we can write :

~( ) ~( ). . ~( ) R f K f H f s s 0 0 0 = ε σ π (8)
where ~( ) H f s 0 is the thermal transfer function defined above.

Two measurements of ~( ) R f are carried out for several frequencies : one on the sample and the other one on a thin reference. For each modulation frequency f i we obtain :

~( )

~( ) . ~( ) R f R f H f s i r i s r s i = ε ε (9) or ~~( ) / R E. H f s r s i = ( 10
)
where E s r = ε ε . E , a , and Bi can be identified by fitting the measured values with the model [START_REF] Siroux | [END_REF].

It is then possible to reach the sample emissivity through two different approaches. The radiometric approach consists in calculating ε s , via E . The calorimetric one consists in calculating via Bi the hemispherical emissivity noted ε s h of the same surface.

In fact, when the sample is in a vacuum chamber, and surrounded by a black surface dimensions of which are greater than those of the sample, the heat transfer coefficient can be written :

h T s h = 4 3 ε σ (11)
where T is the mean temperature of the sample.

So we can obtain ε s h using the below formula :

ε σ s h Bi k T e = . 4 3 (12)
The radiometric approach requires an emissivity reference ε r . On the other hand, the calorimetric approach requires the knowing of the ratio k/e.

Experimental set-up

The apparatus (Fig. 2) is composed of a vacuum chamber (1), a pumping system (2), sample holders (3), an infrared detector (4) and a data acquisition system (5). A turbomolecular pump system permits to reach pressure near 7x10 -5 mbar.

Each sample ( 6) is stock on a modulated sample holder. The latter is composed of a thermoelectric device (3b) squeezed between two metallic plates (3a and 3c). A thermocouple ( 14) is embedded in the metallic plate which holds the sample. Since this plate is thermally thin, readings are assumed to be representative of the front side sample temperature. A cooling water system sinks the heat created by the thermoelectric devices Joule's effect.

The infrared detector is a thermopile with a KRS5 window. The detector is motorized (12) and can be successively placed in front of the different samples.

A PC data based acquisition system controls the excitation voltage of the thermoelectric device, measures sample temperature T(t) as well as the infrared detection system output U(t). The data are first 100 Hz sampled, amplified and analogically low pass filtered (fc = 5 Hz), then they are numerically filtered using a Hanning low pass filter (fc = 0.2 Hz) and 0.4 Hz sampled again.

The excitation voltage is in the form :

V t V V f t mean i i i ( ) .cos( ) = + ∑ - = 2 2 1 0 1 5 π (13)
where V mean is the mean excitation voltage, f 0 is equal to 0.78 mHz. V i is chosen in order that the ratio signal/noise is the same over all the spectrum. The number and the amplitude of each modulation component is chosen regarding the signal/noise ratio, the total peak to peak voltage and the available power. The amplitude of the front side sample temperature is around 12°C peak to peak.

Figure 2.

Experimental set-up : 1-vacuum chamber, 2-vacuum pump system, 3-sample holder (3a, 3c-metallic plate, 3b-thermoelectric device), 4-infrared detector, 5-data acquisition system, 6-sample, 7-mounting sample holder, 8-cooling serpentine, 9-cooling tube, 10-detector holder, 11-sectional iron, 12-motor, 13-computer, 14-thermocouple. For the above configuration (Fig. 2), for a sample the emissivity of which is equal to one, and for a mean temperature equal to 300 K, the minimum theoretical detectable signal amplitude of an infrared detector depends on the detector detectivity and the length of record. It can be shown that (Matteï et al., 1996) for the used thermopile and for L= 85 minutes, this value is around 10 -6 K. However, other noises than detector noise influence to this value, the main one being the D/A board noise which is around 0.05 % of the full scale.

The statistical error of each component of ~( ) U f and ( )

T f 4
is on the order of 0.1 % for the modulus and 0.05 degrees for the phase. These accuracies have been obtained by making measurements of repeatedly and calculating the standard deviation of the resulting distribution of modulus and phase.

If the sample thickness is not very smaller than the lateral dimensions, then the use of a one dimensional model is questionable. In fact, lateral radiant losses influence on back side temperature and identified parameters are then false.

In order to evaluate this influence, a three dimensional numerical model has been developed using finite differences. According to the heat balance equation, the complex temperature Ti of each elementary volume is linked with the complex temperature of the six contiguous elementary volumes :

~~~T T R T jC R T i j j j j i j j j 1 6 1 6 1 (14)
Where R j is a thermal resistance and C i a thermal capacity. Since the problem is symmetric, temperatures are computed for only one quarter of the sample (4608 elements have been used), using an iterative method. The infrared signal of a detector facing the surface sample is then calculated and the Biot number and the diffusivity are identified using a one dimensional analytic model. These values are compared with those obtained with a model without lateral losses, the founded differences are 10 % for the Biot number and 0.1% for the diffusivity. So we can note that the identified value of the Biot number is very dependent on the lateral radiant losses. In order to reduce these losses, reflective surfaces have been closely placed in front of the four sides of the sample.

Results and Discussion

The parameters are identified using the minimization between measured data and calculated ones. In order to take into account the modulus and the phase of data the chosen estimator is the sum of two estimators :

C C C = + ρ ϕ (15) with C i i i k i i k ( ) ( ) * * * 2 1 2 1 , C i i i k i i k ( ) ( ) * * * 2 1 2 1
where and are respectively the modulus and the phase of the data, k is the number of data, * is relating to the measured data, ^ is relating to the calculated data and * is defined by : ϕ ϕ

* * = = ∑ 1 1 k i i k .
Measurements have been carried out on a PVC sample (44mm x 44mm x 5mm) coated with a paint. The diffusivity and the emissivity of the coating are known. The same paint coated on the sample holder is used as the reference.

Figure 4 Identified parameters are : a = 0.113 mm 2 s -1 , Bi = 0.160, E = 1.006.

In order to verify the calorimetric emissivity ε s h we have carried out directional emissivity measurements. Assuming that the normal emissivity is 0.97, calculating calorimetric emissivity using directional emissivity integration yields to a close result ε s h = 0.88. The above results show that it is possible to access to emissivity of thick materials at low temperature via radiometric and calorimetric approaches. This method must be applied now to different sorts of coatings and materials.

Figure 1 .

 1 Figure 1. The conduction model.

Figure 3 .

 3 Figure 3. Amplitude spectra for the temperature and for the infrared signal

Figure 4 .

 4 Figure 4. Modulus and phase of measured values ~( ) / R f s r


O Bokanowski

E Cristiani
email: emiliano.cristiani@gmail.com

H Zidani

An efficient data structure and accurate scheme to solve front propagation problems

Keywords: AMS Classification: 65M06, 49L99 Ultra Bee scheme, narrow band method, sparse matrices, data storage, Hamilton-Jacobi-Bellman equations, front propagation

In this paper, we are interested in some front propagation problems coming from control problems in d-dimensional spaces, with d ≥ 2. As opposed to the usual level set method, we localize the front as a discontinuity of a characteristic function. The evolution of the front is computed by solving an Hamilton-Jacobi-Bellman equation with discontinuous data, discretized by means of the antidissipative Ultra Bee scheme.

We develop an efficient dynamic storage technique suitable for handling front evolutions in large dimension. Then we propose a fast algorithm, showing its relevance on several challenging tests in dimension d = 2, 3, 4. We also compare our method with the techniques usually used in level set methods. Our approach leads to a computational cost as well as a memory allocation scaling as O(N nb ) in most situations, where N nb is the number of grid nodes around the front.

Introduction

This paper is devoted to the study of an efficient numerical technique for coding front propagation in high dimension coming from some optimal control problems.

The idea of level set formulation to propagate curves and surfaces has been introduced in [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] by Osher and Sethian. The advantages of this approach are well known by now. It treats self-intersections, topological changes, kinks, and it is easily extended to capture hypersurfaces in R d , d ≥ 1.

Let us consider a closed moving front Γ t , t ≥ 0, and let Ω t be the closed region that Γ t encloses (i.e, here Ω t is closed and such that Γ t = ∂Ω t ). The level set method aims to find a function u(t, x) such that at each time t the front Γ t can be provided by the 0-level set of the function u(t, •), and we have:

u(t, x) = 0 ⇔ x ∈ Γ t u(t, x) < 0 ⇔ x ∈ Int(Ω t ) u(t, x) > 0 ⇔ x ∈ Ω c t .
1 Moreover, the function x → u(t, x) is uniformly continuous and it increases as the distance between x and Γ t increases. It is known [START_REF] Sethian | Level set methods and Fast Marching methods[END_REF] that u can be obtained by solving a time-dependent Hamilton-Jacobi (HJ) equation with continuous initial datum ϕ vanishing only on the initial front Γ 0 , ϕ < 0 in Int(Ω 0 ), ϕ > 0 in Ω c 0 , and ϕ is a strictly non-decreasing function of the distance to Γ 0 . Many numerical studies have been carried out to construct stable, accurate and efficient methods. See [START_REF] Osher | High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF][START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Shu | High order ENO and WENO schemes for computational fluid dynamics[END_REF][START_REF] Jiang | Weighted ENO schemes for Hamilton-Jacobi equations[END_REF] for a description of such methods on regular grids, and [START_REF] Abgrall | Numerical discretization of first-order Hamilton-Jacobi equation on triangular meshes[END_REF][START_REF] Abgrall | High order numerical discretization for Hamilton-Jacobi equations on triangular meshes[END_REF] on triangular meshes. It is known that in general, once initialized ϕ as the signed distance function to Γ 0 , it is impossible to maintain the level set function as the signed distance function to the moving interface in the advection step. Flat and/or steep regions develop as the interface moves rendering the computation and the localization of the front inaccurate. For this reason, it is necessary to reset the level set function, at regular time intervals, to be the signed distance function to the interface. This step of reinitialization was investigated in [START_REF] Sussman | A level set approach for computing solutions to incompressible 2-phase flow[END_REF][START_REF] Peng | A PDE-based fast local level set method[END_REF] and [START_REF] Hartmann | Differential equation based constrained reinitialization for level set methods[END_REF].

As pointed out in [START_REF] Adalsteinsson | A fast level set method for propagating interfaces[END_REF], an important drawback of the level set approach "stems from the expense; by embedding the interface in R d as the level set of a d + 1-dimensional function, considerable computational labor is required per time step". To overcome this difficulty, Adalsteinsson and Sethian [START_REF] Adalsteinsson | A fast level set method for propagating interfaces[END_REF] suggest a localization of the level set method. This method allows to compute the evolution of the level set function only in a neighborhood around the front. Another fast local level set method has been also proposed by Peng et al. in [START_REF] Peng | A PDE-based fast local level set method[END_REF]. However, in [START_REF] Adalsteinsson | A fast level set method for propagating interfaces[END_REF] as well in [START_REF] Peng | A PDE-based fast local level set method[END_REF] a full d-dimensional matrix is stored, which limits the methods.

In our work, we use a different approach (see Fig. 1). We seek a discontinuous function ϑ(t, x) which takes only values in {-1, 1}

ϑ(t, x) = -1 ⇔ x ∈ Ω t ϑ(t, x) = 1 ⇔ x ∈ Ω c t .
The notion of lower semicontinuous (l.s.c.) viscosity solution [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Frankowska | Lower semicontinuous solutions of Hamilton-Jacobi-Bellman equations[END_REF] leads also to a characterization of ϑ by means of a Hamilton-Jacobi-Bellman (HJB) equation with discontinuous initial datum ϕ given by

ϕ(y) = -1 if y ∈ Ω 0 1 otherwise.
Thus it is natural to concentrate the numerical effort only around the front separating the area of -1 values (Ω t ) from the region of 1 values. For this, we need two ingredients, the first one is a scheme able to compute the front with a good accuracy without diffusion and the second one (as important as the first one) is an efficient way to store and handle only the nodes around the front.

The scheme we consider here is an adaptation of the Ultra Bee scheme used by Desprès and Lagoutière for linear advection and conservation laws [START_REF] Desprès | Un schéma non linéaire anti-dissipatif pour l'équation d'advection linéaire. A non-linear anti-diffusive scheme for the linear advection equation[END_REF][START_REF] Desprès | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF][START_REF] Lagoutière | Modélisation mathématique et résolution numérique de problèmes de fluides compressibles à plusieurs constituants[END_REF][START_REF] Lagoutière | A non-dissipative entropic scheme for convex scalar equations via discontinuous cell-reconstruction[END_REF]. An extension in order to treat HJB equations was proposed in [START_REF] Bokanowski | Anti-diffusive schemes for linear advection and application to Hamilton-Jacobi-Bellman equations[END_REF], numerically showing very good anti-diffusive properties in one and two dimensions (see also [START_REF] Bokanowski | An anti-diffusive scheme for viability problems[END_REF] for computation of the Capture Basin). Convergence results for 1-dimensional HJB equations have been obtained in [START_REF] Bokanowski | Convergence of a non-monotone scheme for Hamilton-Jacobi-Bellman equations with discontinuous initial data[END_REF][START_REF] Bokanowski | Convergence of a non-monotone scheme for Hamilton-Jacobi-Bellman equations with discontinuous initial data[END_REF], where an error bound of order ∆x (the spatial mesh step) in L 1 norm for general bounded l.s.c. initial data is obtained, also proving anti-diffusive properties in particular cases.

In the 2-dimensional case, the Ultra Bee scheme computes values V n i,j which approximate 1 ∆x∆y I ij ϑ(t n , x, y) dxdy, where I i,j is a cell associated to (i, j) and t n is a discrete time. Since ϑ takes values in {-1, 1}, the discrete approximations V n i,j will belong to the interval [-1, 1]. The full matrix (V n ij ) is represented in Fig. 2(left), where values are drawn in white for -1, black for 1, and gray for intermediary values. However, instead of coding the matrix (V n ij ) we can use a sparse matrix where only the values V n ij ∈] -1, 1[ are coded, as well as their first neighboring values, leading to a sparse structure as illustrated in Fig. 2(right). In this paper we use an adapted sparse matrix structure in order to store the discrete front. Note that the idea of using a sparse structure was tested for instance in [START_REF] Megdich | Méthodes anti-dissipatives pour les equations de Hamilton-Jacobi-Bellman[END_REF] for d = 2, 3 (a quadtree structure was also tested in [START_REF] Bokanowski | An adaptative antidissipative method for optimal control problems[END_REF]). However, in these attempts the CPU time is not growing as O(N nb ), where N nb is the number of nodes in the "narrow band" around the front.

In the literature, several techniques for storing sparse matrices have been studied [START_REF] Barrett | Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition[END_REF], [START_REF] Robins | Robs algorithm[END_REF]. These techniques are usually designed to handle standard algebraic computations (for instance, product of matrices). For our purposes, a "good" storage technique should be fast for finding a nonempty value as well as its first neighboring values, and for adding (or removing) nonempty values in the sparse matrix.

Here we shall propose a new storage technique for sparse matrices, and compare it with known storage techniques strategies.

While the worst case complexity of the scheme is no better than dense implementations, in practice it often achieves significant savings O(N nb ), for both memory and CPU time. The advantage of the sparse storage technique we propose is that it can handle larger data sets than other known methods.

The combination of the Ultra Bee scheme and the sparse storage allows us to obtain an efficient algorithm that is able to treat "high-dimensional" problems. We also obtain on some specific examples better results than a traditional level set solver. Finally we mention that the algorithm is quite simple and does not need special features such as reinitialization.

The paper is organized as follows. In Section 2 we recall the HJB equation for ϑ and the Ultra Bee scheme. Section 3 is devoted to the presentation of the sparse dynamic data structure, and to extensive comparisons with other techniques. Section 4 is devoted to numerical examples for the HJB equation. Comparison between the discontinuous approach and the level set approach will be also given.

More applied examples will be treated in a forthcoming work.

Preliminaries

Motivations and setting of the problem

Let Ω 0 be a compact set of R d , and f : R + × R d × R m → R d be a continuous function. We are interested by fast and efficient numerical methods for Hamilton-Jacobi equations of the form:

ϑ t (t, x) + max a∈A(x) {-f (t, x, a) • ∇ϑ(t, x)} = 0, t > 0, x ∈ R d ϑ(0, x) = ϕ(x), x ∈ R d ( 1
)
where ϕ is given by

ϕ(x) = -1 if x ∈ Ω 0 1 otherwise,
and A(x) is a non-empty compact subset of R m , for m ≥ 1. The function ϑ is obviously discontinuous. Nevertheless, the solution of ( 1) is well defined by the notion of lower semicontinuous solution introduced by Barron and Jensen [START_REF] Barron | Semicontinuous viscosity solutions for Hamilton-Jacobi equations with convex Hamiltonians[END_REF] (see also [START_REF] Frankowska | Lower semicontinuous solutions of Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]).

At each time t > 0, the boundary of the set of -1 values of ϑ(t, •) represents the front Γ t which results from the propagation of the initial front Γ 0 := ∂Ω 0 .

Equation (1) models several problems. Let us see some examples.

Advection equation. The simplest situation consists in taking

A(x) = {a}, and f (t, x, a) = f (t, x).

Then equation (1) can simply be rewritten as:

ϑ t (t, x) + f (t, x) • ∇ϑ(t, x) = 0, t > 0, x ∈ R d , ϑ(0, x) = ϕ(x), x ∈ R d .
Here the region Ω 0 is advected with the velocity f . where F is a positive function, then we get the Eikonal equation

ϑ t (t, x) + F (x) ∇ϑ(t, x) = 0 t > 0, x ∈ R d ϑ(0, x) = ϕ(x) x ∈ R d .
The framework of (1) includes also the case of anisotropic propagation. However, it does not include propagation problems with motion by mean curvature.

Target problem: Capture Basin (or backward reachable set). Now we assume that C := Ω 0 is a target. We consider a dynamical system described by the following differential equation: ξ(t) = f (ξ(t), a(t)), for a.e. t ≥ 0,

ξ(0) = x, (2)
where a ∈ L ∞ ([0, ∞[; A) and A is a compact subset of R m (with m ≥ 1). The Capture Basin (or backward reachable set) associated to C is defined as the set of all the initial conditions x ∈ R d from which it is possible to find an admissible trajectory ξ x , solution of (2), reaching the target C before time t ≥ 0:

Capt t (C) := {x ∈ R d : ∃(ξ x , a) satisfying (2) and ∃ 0 ≤ τ ≤ t s.t. ξ x (τ ) ∈ C}.
We define the reachability function ϑ(t, x) by

ϑ(t, x) := -1 if x ∈ Capt t (C) 1 otherwise (3)
Under usual assumptions1 on f , we can characterize ϑ as the unique lower semicontinuous solution of the following HJB equation (by adapting the arguments of [START_REF] Mitchell | A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games[END_REF], see also [START_REF] Megdich | Méthodes anti-dissipatives pour les equations de Hamilton-Jacobi-Bellman[END_REF])

   ϑ t (t, x) + max 0, max a∈A {-f (x, a) • ∇ϑ(t, x)} = 0, t > 0, x ∈ R d ϑ(0, x) = ϕ(x), x ∈ R d (4)
where

ϕ(x) := -1 if x ∈ C, 1 otherwise.
Moreover, from ϑ it is possible to recover the minimal time function T defined by

T (x) := min{t ≥ 0 : ∃ a ∈ L ∞ ([0, t]; A) s.t. ξ x (t) ∈ C}
by means of the following expression

T (x) = min{t ≥ 0, ϑ(t, x) = -1},
with the convention that T (x) = +∞, whenever {t ≥ 0, ϑ(t, x) = -1} = ∅. This result is straightforward and was already stated in the continuous level set setting [START_REF] Osher | A level set formulation for the solution of the Hamilton-Jacobi equations[END_REF].

Rendez-Vous problem. Suppose we want to find out if we can reach the target exactly at time t and not only "before time t". Then instead of (3), we can consider the following definition:

ϑ(t, x) := -1 if ∃ a ∈ L ∞ ([0, t], A) s.t. ξ x (t) ∈ C 1 otherwise. ( 5
)
It is proved in [START_REF] Frankowska | Lower semicontinuous solutions of Hamilton-Jacobi-Bellman equations[END_REF] (see also [START_REF] Barron | Semicontinuous viscosity solutions for Hamilton-Jacobi equations with convex Hamiltonians[END_REF]) that ϑ is the solution of the following HJB equation:

ϑ t (t, x) + max a∈A {-f (x, a) • ∇ϑ(t, x)} = 0, t > 0, x ∈ R d ϑ(0, x) = ϕ(x) x ∈ R d , (6)
which is still is a particular case of (1).

Remark 2.1 For time optimal control problems , in order to reconstruct the optimal trajectories (see [4, Appendix by Falcone]) it is necessary to know the minimal time function T on all over the domain, and not only around the front at the final time. In practice we can compute the solution ϑ together with the minimum time function T : the values of T need only to be saved on the hard disk during computation, before advancing the narrow band. These values are not needed to compute ϑ, and their storage on the hard disk do not reduce the maximal admissible size of the problem which is limited by the RAM's size. (On the contrary, an approach which needs to store the full matrix to make computation will be limited by the RAM memory.) On going work based on the present approach (see [START_REF] Bokanowski | Hamilton-Jacobi-Bellman approach for the climbing problem for heavy launchers[END_REF]) focuses on minimal time function computation and optimal trajectory reconstruction for an optimal control problem including state constraints.

Ultra Bee scheme

We recall here the Ultra Bee (UB) scheme for solving the HJB equation ( 1). This scheme was first studied for advection equations with constant velocity [START_REF] Desprès | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF] (in this context, the scheme is exact). A generalization of the scheme to advection equations with changing-sign velocity is suggested in [START_REF] Bokanowski | Anti-diffusive schemes for linear advection and application to Hamilton-Jacobi-Bellman equations[END_REF], where the properties and the convergence result are proved in dimension 1.

The adaptation of the scheme to solve HJB equations is done in [START_REF] Bokanowski | An anti-diffusive scheme for viability problems[END_REF] (some convergence results are proved in [START_REF] Bokanowski | Convergence of a non-monotone scheme for Hamilton-Jacobi-Bellman equations with discontinuous initial data[END_REF]). Now we present directly the algorithm in dimension 2.

Let ∆t > 0 be a constant time step, and t n := n∆t for n ≥ 0. Let ∆x > 0 be a step size of a spatial grid, and let ξ i,j := (x i , y j ) := (i ∆x, j ∆x) denote a uniform mesh, with i, j ∈ Z. Let us also define

x i+ 1 2 := (i + 1 2 ) ∆x, y j+ 1 2 := (j + 1 2
) ∆x and

I ij :=]x i-1 2 , x i+ 1 2 [×]y j-1 2 , y j+ 1 2 [.
The Ultra Bee scheme aims at computing a numerical approximation of the averages V n i,j := 1 ∆x 2 I ij ϑ(t n , ξ) dξ, for i, j ∈ Z. Since the function ϑ(t n , •) takes only values in {-1, 1}, their averages V n i,j contain the information of the discontinuities localization. The UB-HJB scheme, for the discretization of (1), takes the following form (when there is no ambiguity, we shall omit the time dependence in f and denote

f (ξ ij , a) instead of f (t n , ξ ij , a)) : V n+1 i,j -V n i,j ∆t + max a∈A -f (ξ i,j , a)[D UB a V n ] ij = 0, (7a)
with the initialization

V 0 i,j := 1 ∆x 2 I ij ϕ(ξ)dξ, (7b)
where [D UB a V n ] will play the role of a consistent approximation of the term ∇ϑ(t n , •). To define precisely this approximation, let us first introduce, for j ∈ Z, the fluxes W R j-1 2 and W L j+ 1 2 for a one-dimensional vector (W j ) j∈Z and for real numbers (µ j ) j∈Z ∈ [-1, 1], as follows:

• If µ j ≥ 0, set W L j+ 1 2 :=    min(max(W j+1 , b + j ), B + j, ) if µ j > 0 W j+1 if µ j = 0 and W j = V n j-1 W j if µ j = 0 and W j = V n j-1 , • If µ j ≤ 0, set W R j-1/2 :=    min(max(W j-1 , b - j ), B - j ) if µ j < 0 W j-1 if µ j = 0 and W j = V n j+1 W j if µ j = 0 and W j = V n j+1 , where b + j , b - j , B + j and B - j are defined in (9a)-(9b)). • If µ j ≥ 0 and µ j+1 > 0, set W R j+ 1 2 := W L j+ 1 2 . • If µ j+1 ≤ 0 and µ j < 0, set W L j+ 1 2 := W R j+ 1 2 .
• If µ j < 0 and µ j+1 > 0, then set

W R j+ 1 2 := W j+1 if W j+1 = W j+2 W j otherwise and W L j+ 1 2 := W j if W j = W j-1 W j+1 otherwise. ( 8
) With if µ j > 0,      b + j := max(W j , W j-1 ) + 1 µ j (W j -max(W j , W j-1 )), B + j := min(W j , W j-1 ) + 1 µ j (W j -min(W j , W j-1 )), (9a)
if µ j < 0,      b - j := max(W j , W j+1 ) + 1 |µ j | (W j -max(W j , W j+1 )), B - j := min(W j , W j+1 ) + 1 |µ j | (W j -min(W j , W j+1 )). ( 9b
)
We also use in the sequel the notations F R and F L defined by:

F L (W, µ) j+ 1 2 := W L j+ 1 2 , and F R (W, µ) j-1 2 := W R j-1 2 .
For every (i, j) ∈ Z × Z and every a ∈ A(ξ i,j ), we define

ν 1 i,j (a) := ∆t ∆x f 1 (ξ i,j , a), ν 2 i,j (a) := ∆t ∆x f 2 (ξ i,j , a),
as the local CFL number. In the following we assume that the mesh sizes satisfy the following condition:

CFL := max ∆t ∆x |f 1 (ξ i,j , a)| , ∆t ∆x |f 2 (ξ i,j , a)| ≤ 1. ( 10
)
Now, we define the UB-HJB scheme for (1) as follows (see [START_REF] Bokanowski | Anti-diffusive schemes for linear advection and application to Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Bokanowski | Convergence of a non-monotone scheme for Hamilton-Jacobi-Bellman equations with discontinuous initial data[END_REF]).

Algorithm UB-HJB :

Initialization: We compute the initial averages (V 0 ij ) i,j∈Z as in (7b).

Loop: For n ≥ 0, For a ∈ A, for j ∈ Z, evolve in the x 1 -direction

V n,1 i,j (a) := V n i,j - ∆t ∆x f 1 (ξ i,j , a) F L (V n •,j , ν 1 •,j (a)) i+ 1 2 -F R (V n •,j , ν 1 •,j (a)) i-1 2 , ∀i ∈ Z, where V n •,j = (V n i,j ) i∈Z . Then, for i ∈ Z, evolve in the x 2 -direction V n+1 i,j (a) := V n,1 i,j - ∆t ∆x f 2 (ξ i,j , a) F L (V n,1 i,• , ν 2 i,• (a)) j+ 1 2 -F R (V n,1 i,• , ν 2 i,• (a)) j-1 2 , ∀j ∈ Z, where V n,1 i,• = (V n,1 i,j ) j∈Z . Set V n+1 i,j := min a∈A(ξ i,j ) V n+1 i,j (a) .
Remark 2.2 A general version of the Ultra Bee scheme, for 1-dimensional problems, is given in [START_REF] Bokanowski | Convergence of a non-monotone scheme for Hamilton-Jacobi-Bellman equations with discontinuous initial data[END_REF], for any l.s.c. initial condition in L 1 loc (R). Here, the algorithm is specified to the case of an initial condition taking values only in {-1, 1}.

In [START_REF] Lagoutière | Modélisation mathématique et résolution numérique de problèmes de fluides compressibles à plusieurs constituants[END_REF], the author proved the very interesting property that the Ultra Bee scheme advects exactly a particular class of step functions, in the case of constant advection. We refer to [START_REF] Lagoutière | Modélisation mathématique et résolution numérique de problèmes de fluides compressibles à plusieurs constituants[END_REF][START_REF] Desprès | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF] for other properties.

Storage data structure

In this section we introduce the new data structure we will use with the Ultra Bee scheme to solve numerically equation (1) in high dimension. The goal is to construct a data structure which allows to save memory and to keep acceptable CPU times.

All the numerical results have been obtained on a 2×2GHz Intel processor, 4GB RAM. Codes are written in C++.

Some known storage techniques

First, we recall some classical storage methods. The word "full" will refer to the methods in which every element is stored (using a full matrix), while in the other methods sparse storage techniques are used.

Full Storage (FS): A full d-dimensional matrix is stored, it contains all the values of the function at each grid node. The computation is performed on all over the grid, at each time step. Searching for neighbors is very fast because we have a direct access to each element of the matrix. On the other hand the computation is slow, because at each time step we compute on the full matrix. This is the simplest method to be implemented because there is no need to track the front during its evolution.

Full & List Storage (FLS): this is the classical Narrow Band method (see [START_REF] Adalsteinsson | A fast level set method for propagating interfaces[END_REF]). A full matrix is stored as in FS, but it is also stored a dynamic linked list which contains the indexes of the nodes around the front, that is the zone where we want to compute the solution at each time step. The list must be updated in such a way it follows the front during its evolution. It is not needed to keep the list sorted in any way. Searching for neighbors is done using the full matrix.

List Storage (LS): A dynamic linked list containing the indexes and the values of the nodes around the front is stored, ordered for example in a row-wise fashion with respect to the full matrix (not stored here). Each element of the list has a pointer to the next and to the previous element, so the search for neighbors is done going forward and backward in the list. Moreover, in dimension 2, the indexes (i, j) of a node are "compressed" in only one index k to save memory defining k = iN + j where N is the number of columns. This procedure to "compress" the indexes is easily generalizable to any dimension, provided N and d are not too large (otherwise the integer runs out of bits). This is the best method for memory allocation (see [START_REF] Robins | Robs algorithm[END_REF] for a detailed presentation and comparisons with other methods).

Compressed Row Storage (CRS): This is probably the most common method to store sparse matrices. We use here a slightly different version of the CRS method, substituting static vectors by linked lists to allow resize when needed. As in LS, the full matrix is not stored. See [START_REF] Lin | Efficient data compression methods for multidimensional sparse array operations based on the EKMR scheme[END_REF] for an explanation of the method in dimension d ≥ 2 and an improvement of the basic method for matrix-matrix addition and multiplication.

The Sparse Semi-Dynamic data structure

Now, we introduce what we will call the Sparse Semi-Dynamic (SSD) structure, beginning from the case d = 2. Let us consider without loss of generality a square domain in which the computation is performed. Let M be the N × N matrix which corresponds to the domain. We store a (static) vector p = (p 1 , . . . p N ) of pointers such that every pointer p i is the beginning of a linked list which corresponds to the i-th line of the matrix (not stored), see Fig. 3. Every list is made by elements which contain only the index j and the value of the nodes we want to store in the structure (that is the nodes around the front). The elements are ordered by increasing j's. Note that if a line does not contain any node around the front the corresponding pointer points to NULL. In this way we store only nodes we are interested in and, at the same time, we have a direct access to every line of the matrix, so that we can quickly and easily search for neighbors of a given node. For d = 2 this data structure corresponds to a dynamic version of the CRS method.

In the d-dimensional version of the algorithm the full matrix (not stored) has N d elements. Then we consider a N d-1 -dimensional matrix of pointers. As before, every pointer is the beginning of a list which runs parallel to the last dimension of the matrix, containing all the nodes around the front (see Fig. 4). Note that in the SSD method every element of the lists has only the fields lastindex and value. The first d -1 indexes can be recovered by the pointer we are currently using to run throughout the structure.

For d ≥ 3 our method differs from CRS. To understand the difference, let us focus on the case d = 3. The proposed method, for any couple of indexes (i, j), stores a list corresponding to the Remark 3.1 We can also consider a slightly different storage method, which is a variation of the SSD and that allows to save more memory allocation. The basic idea is that, in the 2D structure of Fig. 3, it is not necessary to store pointers which point to NULL. In this case the vector containing the pointers is substituted by a linked list similar to those containing the values of the nodes. This modification yields to a more complex numerical code. We have verified experimentally that the CPU time is smaller for d = 2, but not always for d > 2.

Comparisons

We now compare the previous storing data structures. We focus on searching neighbors, and on updating the data structures (inserting/removing elements).

Let Γ t be the front at a given time t. We define F + (t) (resp. F -(t)) as the set of nodes which are external (resp. internal) to Γ t with at least one neighbor internal (resp. external) to Γ t .

Here by means of "neighbors" we intend the adjacent nodes in the coordinated directions (i.e. we do not consider the adjacent nodes along the diagonal directions). In all cases but FS and FLS the data structure contains only the nodes in F + ∪ F -. We set the value of the node in F ± equal to ±1. All other values are set to ND (Not Declared).

We will consider two different configurations for the initial front Γ 0 (see Fig. 5):

• Example 1 : Γ 0 is composed by 6 circles (spheres for d = 3, hyperspheres for d = 4) one inside the other. This case is interesting because of the large number of nodes in F + ∪ F -.

• Example 2 : Γ 0 is a rectangle (cylinder for d = 3, hypercylinder for d = 4) with the longest dimension aligned to the d-th one.

All the tests are performed on uniform Cartesian grids on the box [-2, 2] d . We consider the same number N of points in each axis (N d is the total number of grid nodes). We denote by N nb the number of nodes in F + ∪ F -(nodes in the narrow band).

Comparisons for searching neighbors. In this test we measure the CPU time (in seconds) needed to search the value of the 2d neighbors (with d = 2, 3, 4) of all the nodes in F + ∪ F -. The position of the front is fixed (so the data structure is not updated). Results are reported in Tables 1 and2. We can see that FLS is the fastest method because computation is restricted to a narrow band around the front and the research of neighbors is done by using the full matrix. However FLS cannot be used when the number of grid nodes is too large.

In dimension 2, the SSD method coincides with CRS, and their performances are very similar to that of LS. However, as the dimension increases, the difference between these methods becomes larger and larger and the SSD seems to be the only sparse method able to find neighbors in reasonable time.

In Example 1 (see Table 1), we remark that the cost for searching neighbors with the SSD method is of the expected order O(N nb ) = O(N d-1 ). For Example 2, we have numerically observed that the scaling factor of the SSD method remains better then the one of LS or CRS. However, this is the worst situation for SSD in the sense that it does not behave as O(N d-1 ).

On the other hand, if it is known a priori that the front has a dimension predominant over the others, by swapping two coordinates we could easily avoid the worse case for SSD. and remove all the nodes in F -. This simulates an update step for an expanding front. FS is not considered here because it does not need an update procedure. In FLS we do not keep the list sorted, we just add the new nodes at the beginning of the list. The full matrix is used to check if a node to be inserted is already in the list. We perform the tests on the two examples as before, in the case d = 3. Results are summarized in Table 3.

The FLS is clearly the fastest method but it needs much more memory than the sparse methods. SSD greatly overcomes the LS and CRS methods in any test, being only two times slower than FLS in normal situations (Example 1). We shall see in Section 3.5 a concrete example where SSD compares to FLS.

Comments on CPU time and memory allocation

In the SSD structure, although the search for a neighboring value can have a cost up to O(N ) (in the worst case), as well as for updating (adding or removing a node at a given location), we observe in general an order O(1) for a "standard" front propagation problem. Therefore for a narrow band of size N nb = O(N d-1 ), we obtain the CPU time of order O(N d-1 ) in most cases. This will be also exemplified in the following tests (see also next section).

Regarding the size of the front, FS and FLS require a memory allocation of order O(N d ), while LS, CRS and SSD are all of order O(N nb ), which is in general O(N d-1 ).

More precisely, FS and FLS require to store, during all the computation process, a N d -matrix for the values (V n i ) which are real variables (often double precision), while SSD requires only to store a N d-1 -matrix of pointers plus N nb pointers, N nb integers and N nb reals in the dynamic part of the structure.

We also see the memory limitation of the full methods: for d = 3, we could not go further N = 500 approximately (resp. N = 110 for d = 4). For instance, the memory allocated by the different methods, in Example 1 for d=3 and N =400, is: FLS = 565 MB, LS = 54 MB, CRS = 54 MB and SSD = 59 MB.

The Sparse Semi-Dynamic Storage for HJB equations

When using a data structure which does not contain all the nodes of the grid we have to be sure that either the data structure contains all the nodes used for computation or we can recover the values of the nodes even if they are not stored in the data structure. We know that the function ϑ takes value in the set {-1, 1}, so the numerical solution V takes value in the set [-1, 1] because it represents an average of the function ϑ. So, at every time step n, the front is localized in the set of nodes {I : V n I ∈ (-1, 1)} and around this region there is a 1-node thick band where V n = -1 (internal to the front) and a 1-node thick band where V n = 1 (external to the front), see Fig. 2. Only these nodes are stored in the data structure, but they are not enough to perform computation with the UB-HJB scheme. The not stored values can be recovered by the existing ones by the following strategy: if a node is both a neighbor of a node with value -1 (resp., 1) and it is not stored in the data structure, then its value is -1 (resp., 1). Let D n be the set of nodes stored in the data structure at time step t n . Let I be a generic node of the grid and N eigh(I) be the set of the 3 d -1 neighbors of the node I (from now on we have to include also the diagonal directions as required in the UB-HJB scheme). At each time step, D n must first be updated in order to follow the evolution of the front. Then we set

D n+ 1 2 := D n ∪ N eigh(D n ).
In order to maintain the data structure as slim as possible, after each computation (at every time step) we remove from the data structure the set of nodes whose value is either -1 or 1 and it is equal to the values of the neighbors:

R n := {I : V n I = -1 or 1, and

V n I = V n J , ∀J ∈ N eigh(I)}.
Thus,

D n+1 := D n+ 1 2 \ R n .
Now we apply the new data structure to a test problem, solving a simple HJB equation in the form of (4). We set d = 3 and f (x, y, z, a) = (a, 1, 1) T . For testing, we consider N a = 11 discrete control variables a i ∈ [-1, 1] defined by a i = -1 + 2i/(N a -1), i = 0, . . . , N a -1. The initial front (target) is the unit ball. Computation is performed in the box [-2, 2] 3 . The time step ∆t is computed so that the CFL number is 0.9.

The results are summarized in Table 4. We report the CPU time for one time step, without considering the initialization of the data structures nor saving results on disk. |D 0 | is the number of nodes used to locate the initial front. (We have excluded here FS because it is clearly not In this test, the SSD method shows an excellent behavior, with correct scaling, being able to save memory without loosing too much rapidity. The FLS is the fastest method, but it is out of memory as soon as the number of nodes is greater than 500 3 . Also we have numerically observed that as the number of controls N a increases, SSD and FLS have a tendency to have similar CPU times. This is because the cost of managing the sparse structure becomes negligeable with respect to the computational cost for the scheme and which is similar for both approaches.

Numerical examples

In this section we test the UB-HJB scheme with the SSD structure on some numerical examples.

In all tests, we use a variable time step ∆t which is computed at each iteration to fit the CFL condition [START_REF] Bokanowski | An adaptative antidissipative method for optimal control problems[END_REF] with respect to the nodes currently involved in the narrow band. We have chosen CFL=0.9 in all cases.

We shall also compare the results with the level set method based on a very simple RK2-ENO2 solver (without reinitialization), following Osher and Shu [START_REF] Osher | High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF], and which is a priori second order in space and time 2 . Even if this is certainly not the best level set solver we think it is sufficient to illustrate the advantages (or drawbacks) of our UB-HJB discontinuous approach with respect to a continuous level set approach, in particular when the front is not smooth everywhere.

Example 1. Advection of a square. In this first example we solve the following twodimensional advection equation

ϑ t (t, x) + f (x) • ∇ϑ(t, x) = 0, t ∈ [0, T ], x ∈ R 2 , (11)
with a square box initial datum . This is then coupled with a Heun scheme in time (a second order Runge-Kutta scheme).

ϕ U B (x) := -1 if x ∈ [0.5, 1.5] × [-0.
ϑ t (t, x) + ϑ x 1 (t, x) + |ϑ x 2 (t, x)| = 0, t ∈ [0, T ], x = (x 1 , x 2 ) ∈ R 2 . (13a) ϑ(0, x) = ϕ r (x), (13b)
where the initial data is given by

ϕ r (x) := -1 if ||x|| ∞ ≤ r, 1 otherwise (13c)
(here we use the notation

||(x 1 , x 2 )|| ∞ := max(|x 1 |, |x 2 |)).
We consider two types of target:

• r = 0.1: large target case.

• r = 0: thin target case. For numerical purpose, we set r = ∆x (so that there is only one single node that contains a negative value).

Numerically, the equation is discretized on [-1, 2] 2 . Results are given in Tables 7-8 and Figure 6. We notice that, in the large target case, the level set method does not work for N = 51 (i.e. N small), because the set of points such that V n,LS ij < 0 vanishes after a few time steps. The problem is coming from the initial data discretization, which has only one negative value (this is the same when r is only of the order of a few ∆x). This problem would always occurs for thin target problems. Example 3: Two-dimensional deformation of a half plane. We consider a front propagation problem, where the initial front Γ 0 is given by: Γ velocity of the front evolution is given by

0 := {x = (x 1 , x 2 ) ∈ R 2 | x 2 = 0}. The
f (t, x 1 , x 2 ) = -sign( T 2 -t) max(1 -x 2 , 0) -2πx 2 2πx 1 .
The front propagation leads to the following advection equation

ϑ t (t, x) + f (t, x) • ∇ϑ(t, x) = 0 x ∈ R 2 , t ∈ [0, T ], ϑ(0, x) = ϕ(x) (14)
where the function ϕ is continuous in the level set approach and discontinuous in our {-1, 1} approach ϕ LS (x 1 , x 2 ) := min(max(x 2 , -1), 1) in the level set approach

ϕ U B (x 1 , x 2 ) := -1 x 2 ≤ 0 1 otherwise in the {-1, 1} approach. ( 15
)
It is not difficult to prove that the exact solution is given by: ϑ(t, x) = ϕ(R -2πa(x)t 1 x) where R t = cos(t)sin(t) sin(t) cos(t) , a(x) = max(1 -||x|| 2 , 0), and t 1 = min(t, Tt).

Results are given in Figs. 7-8 and Table 9, for T = 6, at two different times: t = 3 (three turns), and t = 6 (return to initial condition after three turns). Exact solution can be obtained by using the method of characteristics. In Figs. 7-8 we show the results for the UB-HJB and level set method, using N = 100 and N = 50 nodes per direction. Anew, the UB-HJB scheme gives good results even on a coarse grid.

Example 4: A three-dimensional rotation problem. Now we consider a three-dimensional advection example:

ϑ t (t, x) + f (x) • ∇ϑ(t, x) = 0, t ∈ [0, T ], x ∈ [-2, 2] 3 , ϑ(0, x) = ϕ(x) (16)
and the corresponding target problem (4), with f (x 1 , x 2 , x 3 ) = (-2πx 2 , 2πx 1 , -1) T where the . For both problems (advection and target), we use only the UB-HJB approximation scheme, see Fig. 9. The computational domain is [-2, 2] 3 . CPU times and errors are summarized in Table 10 (with T = 2 for the advection problem, and T = 0.5 for the target problem). The scaling of the CPU time approximately follows the theoretical scaling of 8 as the mesh size doubles in each direction (a factor 4 comes from the number of cells localizing the surface, and a factor of 2 comes from the time step ∆t, that is divided by 2 because of the CFL condition). The Hausdorff distance decreases well in O(∆x) too.

Here we would have been limited to N = 500 if we had used the FLS approach.

Advection Target problem

(which corresponds to a small square box centered at (0.5, 0, 0.5, 0) T ). For T = 1 the exact solution corresponds to the initial data. Computations are done in the unit box [-1, 1] 4 .

Here we would have been limited to N = 110 (approximately) if we had used the FLS approach. We summarize the numerical results in Table 11. In this example, for N = 200 there is an average of 230, 000 cells in the narrow band (instead of N 4 = 1, 600, 000, 000), 384 time iterations, with an average CPU time of 6.0 seconds per time iteration. Example 6. We consider the HJB equation    ϑ t (t, x) + max 0, max a=1,...,4 (f (x, a) • ∇ϑ(t, x)) = 0, t ∈ [0, T ], x ∈ R 4 , ϑ(0, x) = ϕ(x) [START_REF] Jiang | Weighted ENO schemes for Hamilton-Jacobi equations[END_REF] where the dynamics f (x, a) can have 4 different values: (1, 0, 0, 0) T , (1, 1, 0, 0) T , (1, 0, 1, 0) T or (1, 0, 0, 1) T . The initial data is given by ϕ(x) := -1 if x = (0, 0, 0, 0), 1 otherwise.

The exact solution is given by ϑ(t, x) = -1 on Ω t and ϑ(t, x) = 1 otherwise, where

Ω t := x ∈ R 4 | x 1 , x 2 , x 3 , x 4 ≥ 0, x 2 + x 3 + x 4 ≤ x 1 ≤ t .
The computation is done in [-2, 2] 4 . Results are given in Table 12.

In this example, the estimated CPU time for the computation for one control (one evaluation of the UB scheme for a given control a and for one time step), is about 10% of the CPU time needed to manage the sparse structure for the same time step.

Figure 1 :

 1 Figure 1: Discontinuous approach (left) vs. level set approach (right)

Figure 2 :

 2 Figure 2: Full (left) and sparse (right) matrix coding.

Eikonal equation.

  Let us denote by B(0, 1) the closed unit ball of R d . If we set A(x) = B(0, 1), and f (t, x, a) = F (x)a,

Figure 3 :

 3 Figure 3: a sparse semi-dynamic linked structure storing the nodes

Figure 4 :

 4 Figure 4: SSD structure in three dimensions

10

 10

Figure 5 :

 5 Figure 5: Fronts for example 1 (left) and example 2 (right).

Figure 6 :

 6 Figure 6: (Example 2) Rendez vous problem with N = 101.

6 Figure 7 : 6 Figure 8 : (Example 3 )

 67683 Figure 7: (Example 3) UB-HJB and level set methods at times t = 3 (after three turns) and t = 6 (return to initial data), CFL= 0.9, with N = 100.

Figure 9 :

 9 Figure 9: (Example 4) advection problem (left), target problem (right). T = 2, N = 50.

Table 1 :

 1 CPU times for neighbor searching, Example 1.Comparisons for updating the data structure. Now we compare the CPU time needed to update the data structure. This step consists in inserting and removing nodes, without searching for neighbors. More precisely, we insert all the neighbors of the nodes in F + having value ND

	d	N d	N nb	FS	FLS	order	LS	order CRS order SSD order
	2	400 2	7104	1.2e-3	1.1e-4		2.6e-3	2.8e-3	2.8e-3
	2	800 2	14232	5.0e-3	3.4e-4	3.1 5.3e-3	2.0 5.7e-3	2.0 5.7e-3	2.0
	2 1600 2	28488	2.0e-2	1.3e-3	3.8 1.1e-2	2.1 1.2e-2	2.1 1.2e-2	2.1
	3	100 3	104576	1.3e-2	3.2e-3		1.6e0	2.4e0	2.0e-2
	3	200 3	423704	1.1e-1	2.0e-2	6.2 1.3e1	8.1 1.9e1	7.9 8.3e-2	4.0
	3	400 3	1702192	9.3e-1	9.0e-2	4.5 1.0e2	7.7 1.5e2	7.9 3.4e-1	4.2
	3	800 3	6826816 out of mem. out of mem.	->1e3	->1e3	-1.8e0	4.1
	4	25 4	78704	7.2e-3	2.7e-3		4.2e0	8.1e0	2.1e-2
	4	50 4	674064	1.2e-2	2.6e-2	9.6 1.5e2	35.7 3.1e2	38.2 1.9e-1	9.0
	4	100 4	5549888	2.1e0	3.1e-1	11.9 >1e3	->1e3	-1.4e0	7.1
	4	200 4 45098160 out of mem. out of mem.	->1e3	->1e3	-1.2e1	8.4

Table 2 :

 2 CPU times for neighbor searching, Example 2.

	d	N d	N nb	FLS	order	SSD order
	3 100 3	13560	4.3e-4		6.6e-3
	3 200 3	55808	2.2e-3	5.1 4.5e-2	6.8
	3 400 3	222480	6.6e-3	3.0 3.3e-1	7.3
	3 800 3	894056 out of mem.	-3.0e 0	9.0
	4	25 4	4392	1.4e-4		2.8e-3
	4	50 4	35776	1.0e-3	7.1 3.4e-2 12.1
	4 100 4	297200	1.0e-2	10.0 4.8e-1 14.1
	4 200 4 2450144 out of mem.	-5.9e 0 12.3

Table 3 :

 3 CPU times for updating the data structure, 3-dimensional examples.

	Ex.	N 3	N nb FLS	order LS	order CRS	order SSD	order
	1	50 3	25672 5.3e-3	1.8e-1	2.3e-1	1.1e-2
	1	100 3	104576 2.3e-2	4.3 1.3e0	7.2 2.0e0	8.7 4.9e-2	4.4
	1	200 3	423704 1.0e-1	4.3 1.0e1	7.7 1.7e1	8.5 2.2e-1	4.5
	1	400 3 1702192 4.5e-1	4.5 8.1e1	8.1 1.3e2	7.6 8.9e-1	4.0
	1	800 3 6826816 out of mem.	-6.6e2	8.2 >1e3	-3.8e 0	4.3
	2	50 3	3272 5.0e-4	1.3e-3	3.3e-3	1.3e-3
	2	100 3	13560 2.0e-3	4.0 9.6e-2	7.3 3.7e-2 11.2 6.6e-3	6.2
	2	200 3	55808 8.9e-3	4.4 1.3e-1 13.5 4.9e-1 12.9 3.9e-2	7.1
	2	400 3	222480 3.7e-2	4.1 1.8e0	13.6 5.8e0	12.0 3.0e-1	7.4
	2	800 3	894056 out of mem.	-4.2e1	24.1 7.3e1	12.5 2.9e 0	9.8

Table 4 :

 4 HJB equation, CPU time for one time step (seconds). well as CRS and LS because previous tests showed that they are not better than the SSD method.)

	N 3	|D 0 |	FLS	order SSD order
	50 3	8016	0.07		0.08
	100 3	31416	0.23	3.2	0.37	4.6
	200 3	125968	0.96	4.2	1.45	3.9
	400 3	502368	5.75	6.0	6.52	4.5
	800 3 2011032 out of mem.		27.2	4.1
	1600 3 8043808 out of mem.		113.4	4.1
	efficient, as

Table 6 :

 6 (Example 1) f (x) = (-x 2 , x 1 ): error and total CPU time, T = 2π. CPU order L 1 error Haus. Haus./∆x L 1 error Haus. Haus.

				UB-HJB			Level set
	N 2 /∆x
	50 2	0.04		0.192	0.116	1.45	0.236	0.249	3.11
	100 2 0.16	4.0	0.051	0.056	1.40	0.134	0.151	3.77
	200 2 0.70	4.4	0.049	0.042	2.10	0.038	0.094	4.70
	400 2 2.65	3.8	0.025	0.029	2.90	0.015	0.060	6.00
	Example 2: Two-dimensional Rendez-Vous problem. In this example we consider an
	HJB equation of type (1) with f (x, a) := (1, a) T and a ∈ {-1, 1}, that is

Table 7 :

 7 (Example 2) Rendez-Vous problem, T = 1. Large target (r = 0.1). The level set scheme may loose the front (in this case no error is given).

		UB-HJB	Level set
	N 2 L 1 error Haus. L 1 error Haus.
	51 2	0.178	0.052	-	-
	101 2	0.105	0.022	0.101	0.094
	201 2	0.044	0.011	0.008	0.047
	401 2	0.022	0.006	0.006	0.027

Table 8 :

 8 (Example 2) Rendez-Vous problem, T = 1. Thin target (r = ∆x). In this case, the 0-level set is lost after a few time steps.

						UB-HJB	Level set
				N 2 L 1 error Haus. L 1 error Haus.
				51 2		0.166	0.043	-		-
				101 2		0.080	0.031	-		-
				201 2		0.040	0.016	-		-
				401 2		0.020	0.008	-		-
		EXACT							EXACT
	1	Level Set						1	UB-HJB
	0.8	UB-HJB Initial data						0.8	Initial data
	0.6							0.6
	0.4							0.4
	0.2							0.2
	0							0
	-0.2							-0.2
	-0.4							-0.4
	-0.6							-0.6
	-0.8							-0.8
	-1							-1
		-0.5	0	0.5	1	1.5			-0.5	0	0.5	1	1.5
				Large target						Thin target

Table 9 :

 9 (Example 3) deformation of a half plane at t = 3 (three turns) and t = 6 (return to initial data).

	t = 3

Table 10 :

 10 (Example 4) CPU times and errors for the 3d advection problem (with T = 2) and the corresponding target problem (with T = 0.5), UB-HJB scheme.

		Advection			Target problem
	N 3 CPU	L 1 error Hausdorff	N 3	CPU		L 1 error Hausdorff
	50 3	0.22	4.1e-3	2.6e-1	50 3	0.24		2.3e-1	1.8e-1
	100 3 1.00 4.6	2.2e-3	8.0e-2	100 3 1.15	4.8	1.3e-2	8.0e-2
	200 3 7.19 7.2	6.2e-4	4.4e-2	200 3 8.25	7.2	6.7e-2	4.0e-2
	400 3 64.5 8.9	4.8e-4	2.4e-2	400 3 76.0	9.2	3.5e-2	2.3e-2
	800 3 546.0 8.4	2.6e-4	1.2e-2	800 3 763.0 10.0	1.8e-2	1.3e-2
	initial data is a sphere centered at (-1, 0, 1) with radius r = 0.1:
			ϕ(x) :=	-1 if ||x -(1, 0, 1)|| 2 ≤ r, 1 otherwise
	where ||x||

2 := x 2 1 + x 2 2 + x 2 3

Table 11 :

 11 (Example 5) four-dimensional non constant advection of a square box, T = 1, UB-HJB scheme.

	N 4	CPU order L 1 error Haus. Haus./∆x 1
	25 4	0.89		0.0016	0.140	1.7
	50 4	9.29	10.4	0.0021	0.080	2.0
	100 4	129.1	13.9	0.0013	0.050	2.5
	200 4 2316.2	18.1	0.0007	0.025	2.5

Table 12 :

 12 (Example 6) Results for T = 1, UB-HJB scheme.

	N 4	CPU order Haus. Haus./∆x 1
	25 4	0.28	0.120	0.75
	51 4	1.66	5.9 0.039	0.49
	101 4	11.97	7.2 0.029	0.73
	201 4 175.90	14.7 0.015	0.75

(i) The function f : R d × A -→ R d is continuous; (ii) F (x) := {f (x, a), a ∈ A} is convex compact for all x; (iii) There exists co ≥ 0 s.t. sup a∈A |f (ξ, a)| ≤ co(1 + |ξ|); (iv) For every R ≥ 0, there exists LR ≥ 0, such that ∀y, z ∈ B(0, R), sup a∈A |f (y, a)f (z, a)| ≤ LR|y -z|.

The Hausdorff distance between to sets A, B is defined as dH (A, B) := max(δ(A, B), δ(B, A)) where δ(A, B) := maxa∈A d(a, B).

We consider two possible dynamics: either f (x 1 , x 2 ) = (-2, -1) T with T = 1 (constant advection), or f (x 1 , x 2 ) = (-x 2 , x 1 ) T with T = 2π (rotation). Equation [START_REF] Bokanowski | Convergence of a non-monotone scheme for Hamilton-Jacobi-Bellman equations with discontinuous initial data[END_REF] is discretized on the domain [-2, 2] 2 (with boundary conditions ϑ(t, x) = 1 when needed).

We show in Table 5 and 6 the UB-HJB scheme results, with different discretization numbers (the same number of nodes N is used in each direction). CPU times are given in seconds. The error is computed in L 1 norm between the numerical and exact solution (see below), and the Hausdorff distance 3 is also estimated, between the exact and approximate front.

In the case of the level set method, we choose a regular (Lipschitz continuous) ϕ LS which admits the square border as 0-level set,

In order to use a similar way for computing the L 1 error for both UB-HJB and for the level set approach, we define:

is the numerical value of the UB-HJB scheme, N T is the time step corresponding to the final time T and ξ ij denotes the center of the cell (i, j); we define also

denotes the numerical values computed by the level set method. Table 5: (Example 1) f = (-2, -1): error and total CPU time, T = 1. The L 1 error is 0 for the UB-HJB scheme because of an exact-advection property in this particular example.

UB-HJB

Level set N 2 CPU order L We remark that the L 1 error estimate for UB-HJB may vanish, in the constant advection case, which means here that the computed points with value -1 are all in the expected correct cells. This comes from a known exact-advection property of the Ultra Bee scheme as mentioned in Section 2.2.

We also note that the UB scheme, which is only a first order scheme at best, gives here results comparable to the second order RK2-ENO2 level set method. This is because the level set method is second order only in the regions where the front is regular, and is not expected to be better than first order otherwise (the theoretical bound is O(∆x 1/2 ) for monotone schemes, following [START_REF] Crandall | Two approximations of solutions of Hamilton-Jacobi equations[END_REF]).


Kevin Nadaud
email: kevin.nadaud@univ-tours.fr

Guillaume F Nataf

Nazir Jaber

Micka Bah

Béatrice Negulescu

Pascal Andreazza

Pierre Birnal

Jérôme Wolfman

Subcoercive-field dielectric response of 0.5 (Ba 0.7 Ca 0.3 TiO 3 ) -0.5 (BaZr 0.2 Ti 0.8 O 3 ) thin film: peculiar third harmonic signature of phase transitions and residual ferroelectricity

Keywords: Impedance spectroscopy, hyperbolic analysis, relaxor, domain wall

Sub-coercive field non-linearities in 0.5 (Ba 0.7 Ca 0.3 TiO 3 ) -0.5 (BaZr 0.2 Ti 0.8 O 3 ) (BCTZ 50/50) thin film elaborated using pulsed laser deposition are studied using permittivity and phase angle of the third harmonic measurements as function of the AC measuring field E AC and temperature. The global phase transition temperature T max for which the permittivity is maximum, decreases from 330 K to 260 K when E AC increases. Rayleigh analysis of the AC field dependence of the relative permittivity shows a regular decrease of the domain wall motion contributions as temperature increases up to T max and an even more pronounced decrease above T max . This measurement reveals that the ferroelectric behavior subsists 70 K above the global phase transition. The phase angle of the third harmonic at temperatures below 275 K, is characteristic of a conventional ferroelectric and from 275 K to T max = 330 K of a relaxor. Above T max , the thin film exhibits a peculiar phase angle of the third harmonic, which consists of -180 • → -270 • → 0 • instead of the -180 • → -90 • → 0 • found for relaxor.

In this article, we study the phase transitions in a thin film of BCTZ 50/50 using the relative permittivity at the first harmonic (using hyperbolic law) and the evolution of phase angle of the third harmonic as a function of the AC field for different temperatures. Using hyperbolic analysis, we show that residual ferroelectricity persists up to 70 K above the maximum in permittivity indicating the transition to the cubic phase. Above this temperature, a faster decay of the reversible and irreversible domain wall motion contributions is visible. The crossing of this temperature is also well visible on the phase angle of the third harmonic, which exhibits a totally different evolution with increasing AC field below and above it.

At sub-switching AC fields, for a homogeneous distribution of pinning centers, the relative permittivity can be described using the Rayleigh law: 28,29 ε r = ε r-l + α r E AC (1) Where E AC corresponds to the magnitude of the applied measuring electric field, ε r-l is the lattice contribution to the permittivity, α r is the irreversible contribution from the motion of domain walls (domain wall pinning/unpinning), polar cluster boundaries, or phase boundaries and corresponds to the slope of the permittivity v.s. electric field curve. In such conditions of homogeneous distribution of pinning centers, the polarization versus electric field loop can be described using the following expression: 20,28

P = ε 0 (ε r-l + α r E 0 ) E ± ε 0 α r 2 E 2 0 -E 2 + . . . ( 2
)
The sign + stands for the decreasing and the signfor the increasing part of the AC field. The second term reflects the hysteretic contribution of domain walls to the polarization. This non-linear expression of the polarization gives the following Fourier series decomposition when the applied electric field is E(t) = E 0 sin (ωt): 30

P (t, E 0 ) = ε 0 (ε r-l + α r E 0 ) E 0 sin (ωt) + 1,3,5,... 4ε 0 α r E 2 0 sin πn 2 πn(n 2 -4) cos (nωt) (3)
The irreversible domain wall motion contribution is thus out-of-phase with the measuring electric field in the case of an ideal material. In order to describe a real material, equation ( 2) can contain additional terms, reflecting the degree of randomness of the energy profile, 21 and in that case, harmonics may not be purely out-of-phase. For this reason, the non-linear response of a relaxor or a ferroelectric material can be investigated by extracting the phase of the third-harmonic contribution to the polarization and its evolution with the measuring field amplitude, which gives information on the hysteretic or nonhysteretic character of domain wall motions 13,20,21,31,32 . In a real material, the distribution of pinning centers is not homogeneous and for low AC fields, the relative permittivity is almost constant, corresponding to reversible domain wall contributions, also called domain wall vibration 29,33 . A generalized expression can then be used to describe the permittivity evolution, called the hyperbolic law: [34][35][36]

ε r = ε r-l + ε 2 r-rev + (α r E AC ) 2 (4)
with ε r-rev the reversible domain wall motion contribution, proportional to the domain wall density [37][38][39] . ε r-l , ε r-rev and α r can be obtain by measuring the relative permittivity as a function of the driving field. Their evolution with driving frequency, temperature, DC bias field or previous states, allows understanding dielectric relaxations 40 , residual ferroelectricity 18 , contributions to tunability 41,42 and annealing/cycling effects 32,38 . Polycristalline 380 nm thick BCTZ film has been grown on Pt/TiO 2 /MgO substrate by pulsed laser deposition. Details on growth conditions can be found in the supplementary information. Top circular Au/Ti electrodes (150 µm radius) were deposited through a shadow mask. The metal-insulator-metal topology has been chosen for the simple extraction of the dielectric properties using the parallel plate capacitor formula.

X-ray diffraction (see diffraction pattern in Fig. S1 of supplementary information) shows that the BCTZ film is single-phased and polyoriented. Peak positions could be indexed according to orthorhombic BCTZ 50/50 (pdf #04-022-8189). The film composition has been characterized by Rutherford Back Scattering (RBS, see Figure S2a,b) and was found to be close to the nominal target composition (see table S1) .

The single frequency Vertical Piezoresponse Force Microscopy (VPFM) is used to map the electrical polarization at nanoscale in the BCTZ film. A SCM-PTSI tip with a spring constant of 2.8 N m -1 and a radius of curvature of 15 nm is mounted on the probe holder. An AC amplitude of 0.5 V at a drive frequency of 327 kHz (i.e. the contact resonance frequency) is applied to the BCTZ sample while the tip is grounded. The scan size and the scan rate are set to 500 nm and 0.2 Hz, respectively.

The dielectric characterizations presented in this article have been acquired using a lock-in amplifier (MFLI with MD option, Zurich Instrument) connected to a temperature-controlled probe station (Summit 12000, Cascade Microtech). The AC measuring signal has been generated using the embedded lock-in generator. Its amplitude has been swept from 10 mV rms to 1 V rms at a frequency of 10 kHz. The applied voltage and current through the capacitor are measured simultaneously by the lock-in. The first harmonic of the voltage, in addition to first, second and third harmonics of the current, are demodulated simultaneously. |V | exp (jθ V ) is the phasor representing the applied voltage and |I k | exp (jθ I k ) the phasor representing the k-th harmonic of the current (with j the imaginary unit). The first harmonic of the current and the applied voltage are used to compute the complex impedance: 32,35,43

Z = |V | |I 1 | exp (j (θ V -θ I1 )) (5)
The complex capacitance C * can hence be derived from the complex impedance:

C * = 1 jωZ = |I 1 | ω|V | exp j θ I1 -θ V - π 2 (6)
ω the angular frequency of the measuring voltage. The material relative permittivity ε * r :

ε * r = t Sε 0 C * , ( 7
)
with t the thickness of the film, S the surface of the electrodes and ε 0 the vacuum permittivity. In the present case, the electrodes are sufficiently thick in order to limit the effect of the series resistance on the measured impedance.

In addition to the measurement as a function of the AC amplitude, the relative permittivity has been measured from 31 Hz to 31 kHz with an AC amplitude of 40 mV rms .

The k-th harmonic phase angle extraction process is similar and its value can be obtained using the following expression: 32

δ k = θ I k -kθ V - π 2 (8)
Fig. 1a displays the PFM phase, associated with the domain orientation (up or down). It highlights that the majority of the grains in the thin film form a single domain in the out-of-plane projection of the electrical polarization. Nevertheless, few grains present a polydomain configuration, as the superimposition of the topography one and that of the PFM phase shows (see dark and bright colors in Fig. 1b ). Fig. 1c displays the PFM phase and amplitude signals as a function of DC bias. These piezoresponse data are a signature of ferroelectric behavior. A phase difference of 182°between the up and down domains is observed at the coercive voltage (|V c | ≃ 1 V). At the same time, the PFM amplitude drops as the polarization switching under the tip gradually takes place. Topography measurement using AFM is given in supplementary material.

Fig. 2a shows the relative permittivity as a function of temperature for different frequencies. T max have been determined from a least square fitting of ε ′ r with a parabola 44 . The decrease of the permittivity when the frequency increases is stronger for low temperatures (below 300 K) than than for high temperatures (above 340 K). This results in a shift of the maximum permittivity temperature T max with frequency (Fig. 2b), typical of relaxor-ferroelectrics behavior 12,15 . The relaxor behavior is confirmed by the P (E) loop and the analysis of the phase transition diffuseness, using modified Curie-Weiss analysis, both in supplementary material (Fig. S4 andS5).

Fig. 2c shows the relative permittivity as a function of temperature for different values of the AC measuring field. For an AC field of 0.3 kV cm -1 , T max is 330 K, i.e. slightly lower than that reported for bulk BCTZ ceramics in the literature (363 K 45 , 365 K 46 , 373 K 47 ) but very similar to the one obtained for thin films: 323 K 48,49 . An inflection point is also visible, at 270 K, and may correspond to a structural transition. It is close to the tetragonal-orthorhombic and orthorhombicrhombohedral transition temperatures observed in bulk ceramics 50 and matches with anomalies observed in pyroelectric current 51 and dielectric and elastic properties 45 . Bulk 0.5 (Ba 0.7 Ca 0.3 TiO 3 ) -0.5 (BaZr 0.2 Ti 0.8 O 3 ) (BCTZ 50/50) is cubic (Pm-3m) above 360 K. On cooling, it enters a tetragonal phase (P4mm), that coexist with a rhombohedral phase (R3m) between 310 K and 210 K. Below 210 K, it is rhombohedral 50 .

When the AC measuring field increases, the relative permittivity increases due to irreversible contributions from the motion of domain walls, polar cluster boundaries, or phase boundaries [52][53][54] . One can note that the lower the temperature, the higher the increase of the relative permittivity with AC field. The consequence of this increase is a shift towards lower temperatures of the temperature where the relative permittivity is maximum (T max ) when the AC field increases (Fig. 2d). Thus, in the following, the estimation of the maximum of the permittivity T max indicative of the transition to the cubic phase, is conducted at low AC field, to reduce the influence of domain wall motion.

Fig. 3 shows the relative permittivity variation ∆ε ′ r = ε ′ r -ε ′ r (0) and the phase angle of the third harmonic, as a function of the AC measuring field amplitude and for different temperatures. ∆ε ′ r increases when the AC field increases, due to the irreversible domain wall mo- tion contribution. This increase is stronger at low temperatures (Fig. 3a,3c). The increase in permittivity vs AC field persists above T max (black dashed line), corresponding to residual ferroelectricity, similar to what has been observed in (Ba, Sr)TiO 3 17,18,55 . Fig. 3b shows the phase angle of the third harmonic as a function of the AC measuring field, for different temperatures. At low temperatures, up to 275 K, the evolution of the phase-angle is typical of soft or weakly hard ferroelectrics 20 : for low AC fields, δ 3 ≃ -180 • , then when the AC field increases, δ 3 increases and tends to a limit value around -45 • . This transition indicates a flat and random pinning potential for domain walls that can easily move both in reversible and irreversible ways. At low fields, a value of δ 3 = -180 • corresponds to a non-hysteretic contribution of the domain wall motion to the permittivity 13,30 . This value is found for many ferroelectrics at low AC fields where the reversible domain motion contribution (vibration) dominates 29,32,33 . At higher fields, the phase angle of the third harmonic evolves toward -90 • and corresponds to a purely hysteretic domain wall motion contribution to the permittivity, which results in a linear increase of the permittivity. It is the theoretical Rayleigh behavior. At even higher fields, a saturating-like response is observed with δ 3 around -45 • . The evolution of the permittivity for these temperatures (Fig. 3a) also corresponds to what is obtained for soft or weakly hard ferroelectrics: a very limited range for which the permittivity linearly increases (from 3 kV cm -1 to 8 kV cm -1 ) then, the permittivity follows a sublinear increase. When the temperature increases, the phase evolution with the AC field is similar up to 275 K, the main difference resides in the phase-angle value at large fields (33 kV cm -1 ) which progressively increases.

At 300 K, the phase-angle increases from -180 • to -135 • and then increases quickly to 0 • . This type of phase-angle response is close to the one obtained for relaxors 13,20,21,30 hysteretic contribution, but in that case it corresponds to a saturation of the permittivity. The principal difference with the literature is the presence of a plateau before the increase to 0 • which extends to larger electric fields here (20 kV cm -1 instead of few kV cm -1 ). This may be due to the different sample form (thin film here vs ceramics) [similar to Fig. S4 in 21]. At 325 K, the response is quite similar to the one at 300 K, even if the transition from -180 • to 0 • is sharper.

For temperatures above T max , i.e. above 330 K, the phase-angle response exhibits a peculiar evolution: at low field the value stays almost constant at -180 • and instead of a transition by -90 • to go to 0 • , the phaseangle decreases and pass by -270 • . The crossing of T max and the associated change of regime is well visible on Fig. 3d. This behavior has not been reported in the literature and this abrupt change in the phaseangle can be used to monitor the crossing of T max . A phase-angle of -270 • indicates a pinching of the hysteresis loop 31 . For the pinching, a value close to -240 • is usually observed experimentally because of the competing contribution of the vibration from domain walls (δ 3 = -180 • ) 30 . Pinching of the hysteresis loop can have many origins: (i) fresh state for hard ferroelectrics such as Fe-doped PbZr 0.58 Ti 0.42 O 3 31 , (ii) aged state and oxygen vacancies migration for Cu-doped BaTiO 3 56 , Ce-doped Ba(Ti 0.99 Mn 0.01 )O 3 57 or (iii) electric field-induced transition from the paraelectric to a ferroelectric state for BaTiO 3 [START_REF] Merz | [END_REF] or (1x)(Bi 1/2 Na 1/2 )TiO 3x BaTiO 3 59 . This later could play a role in our measurement since the change occurs near T max and the magnitude of the AC field stays low (and limits aging).

The electric field E t1 needed to induced a paraelectric to a ferroelectric state above T max should increase with the temperature [START_REF] Merz | [END_REF]59 . One can note in our case, the electric field E -270 • AC , for which the phase angle of the third harmonic pass by -270 • , decreases from 18 kV cm -1 to 12 kV cm -1 when the temperature increases from 335 K to 380 K. Thus, the field E -270 • AC may not correspond to E t1 but to the field for which (i) below, the non linearity is governed by reversible domain wall contribution (δ 3 ≃ -180 • ) and (ii) above, the non-linearity is governed by saturation (δ 3 ≃ 0 • ). Around E -270 • AC , the deviation from the theoretical value, δ 3 = -180 • for reversible and δ 3 = 0 • for saturation, is attributed to a pinching of the loop and the measured phase corresponds to the combination of the different contributions. 30 The real part of the permittivity measurement as a function of the AC measuring field has been fitted using (4) and Levenberg-Marquat algorithm and the resulting fits are shown as solid lines (Fig. 3a). Only the data below 7 kV cm -1 have been used for the fit since a deviation from the Rayleigh law occurs at higher fields. The hyperbolic fits values are given in Fig. 4 as a function of the temperature, with associated confidence intervals. The lattice contribution shows a similar evolution with temperature as the total permittivity, consistent with the fact that, at low fields, it represents the main contribution to the permittivity. Its maximum occurs for 334 K, very close to the T max found at low fields (330 K). The slight difference below T max correspond to the reversible domain wall motion contribution.

Fig. 4b shows the domain wall motion contributions to the permittivity. There a no measurements in the literature on thin films of BCTZ. Still it is useful to compare with room temperature measurements on ceramics where α r is around 300 cm kV -1 to 3000 cm kV -160-62 , much larger than the values we find here. This indicates that domain walls are less mobile in thin films which is attributed to clamping of the film 63 .

Below T max , both domain wall motion contributions decrease when the temperature increases, without any change at 270 K where there was an inflection point in the real part of the permittivity. This differs from what has been found for films of (Pb, Sr)TiO 3 36 or 0.5 PbYb 1/2 Nb 1/2 -0.5 PbTiO 3 33 for which the domain wall motion contribution increases when the temperature approaches T max . Above T max , the domain wall motion coefficients still decrease but are not null, corresponding to a residual ferroelectricity 17,18 , which persists here 70 K above T max . The decrease is slightly more pronounced than before: at 330 K α ′ r is 3 times lower than at 260 K whereas at 400 K, α ′ r is 9 times lower than a 330 K. The change of decay rate is also visible on the reversible contribution. According to Boser 37 , the reversible domain wall contribution is proportional to the domain wall density and the distance traveled by the domain wall, this later depending on the potential energy profile of the domain wall 33 . A stronger decrease of the reversible contribution could thus indicate faster changes in the density of domain walls, which would be consistent with the fact that above T max the number of PNRs and their sizes both decrease. 44,64 In the studied sample, the reversible and irreversible contributions do not exhibit a peak, contrary to what has been seen for NaNbO . Using reversible and irreversible domain wall motion contributions, it is possible to determine the threshold field E th = ε ′ r-rev /α ′ r which represents the degree of pinning of the domain walls 34,39,42 . The threshold field remains almost constant when the temperature changes, in- dicating a constant degree of pinning. This suggests that the decrease of non-linear behavior when the temperature increases has for origin the decrease of the domain wall density and not the change of the energy profile.

In summary, in this article, we studied the sub-coercive field non-linearities as a function of the temperature of a BCTZ 50/50 thin film. The relaxor behavior of the film is confirmed by the shift of the maximum permittivity temperature T max with frequency and the slim shape of the P (E) loops. We have shown that T max also depends on the AC measuring field because of the irreversible do-main wall motion contribution. The domain wall contributions is found to regularly decrease with temperature and still contribute to the permittivity above T max . In addition, phase angle of the third harmonic measurements show that the thin film behaves like a conventional ferroelectric below 275 K and as a relaxor from 275 K to T max = 330 K. Above T max , the thin film exhibits a peculiar phase angle of the third harmonic response which consists of -180 • → -270 • → 0 • instead of the -180 • → -90 • → 0 • found for relaxor. The change in phase angle of the third harmonic may be used to monitor the crossing of T max .

Figure 1 .

 1 Figure 1. Phase of the Vertical PFM (VPFM) (a). Superimposition of the phase signal and the topography (b). PFM phase and amplitude as a function of an applied DC voltage VDC (c).

  EAC = 33 kV cm 1 EAC = 27 kV cm 1 EAC = 22 kV cm 1 EAC = 16 kV cm 1 EAC = 11 kV cm 1 EAC = 5.9 kV cm 1 EAC = 0.4 kV cm1

Figure 2 .

 2 Figure2. Relative permittivity as a function of the temperature for different measurement frequency for E AC = 1.2 kV cm -1 (a). Tmax for which the relative permittivity is maximum as a function of the frequency for E AC = 1.2 kV cm -1 (b). Relative permittivity as a function of the temperature for different values of the AC measuring field for f = 10 kHz (c). Tmax for which the relative permittivity is maximum as a function of the AC measuring field for f = 10 kHz (d).

Figure 3 .

 3 Figure 3. Relative permittivity variation (a)-(c) and phase angle of the third harmonic (b)-(d) as a function of the AC measuring field for different values of the temperature. For relative permittivity presented in (a), dots correspond to the experimental data and the solid lines to the hyperbolic fits. The black dashed line corresponds to Tmax (at E AC = 0.3 kV cm -1 ). The red dashed line in (c)-(d) corresponds to the inflection point at T = 270 K. Data presented in (a)-(b) corresponds to slices every 25 K of the data presented in (c)-(d).

1 )Figure 4 .

 14 Figure 4. Relative permittivity for a measuring field of 0.3 kV cm -1 , lattice contribution to the permittivity (a), domain wall motion contributions to the permittivity (b) and threshold field (c) as a function of the temperature.

  3 52 , Pb 0.92 La 0.08 Zr 0.52 Ti 0.48 O 3 44 , Nb-doped PbZrO 3 64 or PbMn 1/3 Nb 2/3 O 3 54

ACKNOWLEDGMENTS

This work has been performed with the means of the CERTeM (microelectronics technological research and development center) of French region Centre Val de Loire. This work was funded through the project MAPS in the program ARD+ CERTeM 5.0 by the Région Centre Val de Loire co-funded by the European Union (ERC, DYNAMHEAT, N°101077402). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

CONFLICT OF INTEREST

The authors declare no competing financial interest.