Olivier Ley

AND Erwin Topp

Miguel Yangari

SOME RESULTS FOR THE LARGE TIME BEHAVIOR OF HAMILTON-JACOBI EQUATIONS WITH CAPUTO TIME DERIVATIVE

Keywords: Mathematics Subject Classification. 35K55, 35R09, 47G20, 35B40, 33B20 Integro-Differential Equations, Nonlinear Partial Differential Equations, Caputo derivative, Regularity, Large Time Behavior. 1

HAL is

Introduction.

In this note we are interested in nonlocal Hamilton-Jacobi equations with the form ∂ α t u + H(x, Du) = 0 in Q := T N × (0, +∞). (1.1) subject to the initial condition u(•, 0) = g in T N , (1.2) for some H ∈ C(T N × R N ) and g ∈ Lip(T N ) given.

The nonlocal nature of the problem is cast by the operator ∂ α t , which denotes the Caputo time derivative of order α ∈ (0, 1), starting at time zero. For φ ∈ C 1 (0, +∞) it is defined as

∂ α t φ(t) = 1 Γ(1 -α) t 0 φ (s) |t -s| α ds, (1.3)
where Γ is the Gamma function that acts as a normalizing constant making ∂ α t become the usual first order derivative when α → 1 -, see [START_REF] Diethelm | The analysis of fractional differential equations[END_REF][START_REF] Gorenflo | Mittag-Leffler functions, related topics and applications[END_REF] and references therein. Following the ideas of [START_REF] Allen | A parabolic problem with a fractional time derivative[END_REF], where under appropriate assumption on the function φ, its Caputo derivative in (1.3) can be equivalently computed as (1.4)

∂ α t φ(t) = cα t -∞ φ(t) -φ(s) |t -s| 1+α ds,
for some cα > 0, and where we have extended φ as φ(t) = φ(0) for t < 0. This operator has a nonlocal degenerate elliptic nature in the sense of Barles and Imbert [6] that allows to conclude comparison principle among viscosity solutions as it is proved in [START_REF] Topp | Existence and uniqueness for parabolic problems with Caputo time derivative[END_REF]. This is a powerful reason to consider Caputo derivative instead of other fractional derivatives, as, for example, Riemann-Liouville derivative defined for adequate functions φ as

∂ α RL φ(t) = c α d dt t 0 φ(z) (t -z) α dz.
Moreover, as it can be seen in [START_REF] Diethelm | The analysis of fractional differential equations[END_REF]Chapters 5,[START_REF] Barles | Second-order elliptic integro-differential equations: viscosity solutions' theory revisited[END_REF], Caputo derivative is more adequate to deal with the most classical notion of initial condition compared with Riemann-Liouville problems, where the initial condition is understood in a generalized sense.

Coming back to (1.1), and more specifically to the Hamiltonian H, we assume throughout this paper H(x, p) is periodic in x and coercive in p.

We focus our attention into Bellman-type Hamiltonians with the classical structure related to optimal control problems with compact control set, that is, satisfying the regularity/growth condition for some c H > 0, and for all x, y ∈ T N , p ∈ R N . We are also interested in the case H has superlinear growth in the gradient, common in control problems with unrestricted control space. We refer to this case through the following assumption: there exists m > 1 (the "gradient growth") and A > 0, c H > 1 such that for all µ ∈ (0, 1), x, y ∈ T N , p, q ∈ R N with |q| ≤ 1. We notice at once that the first condition in (1.6) implies that there exist C > 0 such that (1.7) H(x, p) ≥ C -1 |p| m -C for all x ∈ T N , p ∈ R N , and the same inequality holds with m = 1 if the second condition holds in (1.6) and H is convex.

Our interest is the analysis of the behavior of the solutions to (1.1)-(1.2) and make a contrast with the classical case α = 1, namely, the Hamilton-Jacobi equation

(1.8) ∂ t u + H(x, Du) = 0 in Q.
There is a vast literature regarding problem (1.8)-(1.2). We refer the surveys of Barles and Ishii in [START_REF] Achdou | Hamilton-Jacobi equations: approximations, numerical analysis and applications[END_REF] and references therein for the basics about this problem, like existence, uniqueness and regularity. A natural question that arises is the analysis of the behavior of the solution for long times. There too, there are a lot of references when α = 1, see [START_REF] Fathi | Sur la convergence du semi-groupe de Lax-Oleinik[END_REF][START_REF] Namah | Remarks on the long time behaviour of the solutions of Hamilton-Jacobi equations[END_REF][START_REF] Barles | On the large time behavior of solutions of Hamilton-Jacobi equations[END_REF][START_REF] Fathi | Existence of C 1 critical subsolutions of the Hamilton-Jacobi equation[END_REF][START_REF] Davini | A generalized dynamical approach to the large time behavior of solutions of Hamilton-Jacobi equations[END_REF][START_REF] Ishii | Asymptotic solutions of Hamilton-Jacobi equations for large time and related topics[END_REF][START_REF] Barles | A new PDE approach to the large time asymptotics of solutions of Hamilton-Jacobi equations[END_REF][START_REF] Achdou | Hamilton-Jacobi equations: approximations, numerical analysis and applications[END_REF] and the references therein. Here, we focus on the nowadays classical framework of Namah and Roquejoffre paper [START_REF] Namah | Remarks on the long time behaviour of the solutions of Hamilton-Jacobi equations[END_REF], where the authors address the long time behavior of (1.8) under the assumption that (1.9)

          
H(x, p) = F (x, p) -f (x) for all (x, p), with F , f continuous, F (x, •) is convex for all x ∈ T N , F (x, p) > F (x, 0) = 0 for all x ∈ T N , p ∈ R N \ {0},

f ≥ 0 in T N .
It is possible to generalize slightly the above assumptions but we choose to state them in this form since all the main difficulties are present.

Under these assumptions, (1.1) reads

∂ α t u + F (x, Du) = f (x), (1.10) f ∈ Lip(T N ) and Z := {x ∈ T N : f (x) = min T N f }, (1.11)
the so-called projected Aubry set [START_REF] Fathi | Existence of C 1 critical subsolutions of the Hamilton-Jacobi equation[END_REF], which is a compact subset of T N . It follows from Lions, Papanicolaou and Varadhan [START_REF] Lions | Homogenization of Hamilton-Jacobi equations[END_REF] that the so-called ergodic problem

(1.12) F (x, Dv) = f (x) + c x ∈ T N , has a solution (c, v) ∈ R × W 1,∞ (T N
) and c is unique. Actually, under our assumptions, it is easy to see that c = -min T N f . It is then expected that the long time behavior of the solution u of (1.1)-(1.2) is given by the asymptotic expansion u(x, t) + ct α = v(x) + o(1), o(1) → 0 uniformly as t → +∞, (1.13) where (c, v) is a solution of (1.12).

In the local case α = 1, this asymptotic behavior is established in [START_REF] Namah | Remarks on the long time behaviour of the solutions of Hamilton-Jacobi equations[END_REF]Theorem 1]. The proof relies basically on three steps. The first step is to obtain that the set {u(•, t) + ct, t ≥ 0} is relatively compact in W 1,∞ (T N ). Let us point out that this property is not sufficient and the main difficulty is to prove the full convergence of u(•, t) + ct as t → +∞. The second step is to notice that ∂ t (u(x, t) + ct) ≤ 0 for x ∈ Z, from which one infers that u(•, t) + ct is nonincreasing in time on Z, so it converges uniformly to a Lipschitz continuous function φ on Z as t → +∞. The third step is to take the half-relaxed limits

u(x) = lim sup y→x,s→t, →0 u(y, s ) + c s and u(x) = lim inf y→x,s→t, →0 u(y, s ) + c s ,
which are, respectively, a sub-and a supersolution of (1.12) and then to apply a strong comparison result for (1.12) with the "Dirichlet boundary condition" u = u = φ on Z. It follows u = u in T N , which gives the desired full convergence.

Our goal is to develop a similar procedure for the Caputo fractional case α ∈ (0, 1). The first step holds. Indeed, the elliptic properties shown by the expression (1.4) are not only related to comparison principles but also to regularity in the time variable. Using Ishii-Lions method for nonlocal problems presented in [START_REF] Barles | Hölder continuity of solutions of second-order nonlinear elliptic integro-differential equations[END_REF], we are able to prove that bounded solutions to (1.1) are α-Hölder continuous in time, uniformly in Q, see Theorem 2.1. We recall that such a property does not come from the contraction principle given by comparison arguments as in the classical case, basically because such a property is not known for fractional problems, where the influence of the "memory" put troubles in the analysis of problems shifted in time. Moreover, in the local case (1.8), Lipschitz estimates in time allows to use Rademacher's Theorem to regard u t as an L ∞ function, and extract the boundedness of Du through the coercivity of H. Such a program cannot be carried out directly in (1.1) since α-Hölder functions are not sufficient to make ∂ α t u in L ∞ . Nevertheless, we can get equicontinuity in the space variable by a regularization procedure via inf and sup convolutions, see Theorem 2.2.

Similarly, the third step is identical to the one in the classical case since the ergodic problem (1.12) is the same for all α ∈ (0, 1].

It follows that the limiting step is the second one. We still have ∂ α t (u(x, t) + ct α ) ≤ 0 for x ∈ Z but it is not anymore sufficient to infer neither that u(•, t) + ct α is nonincreasing in time on Z, nor that it converges. As we show in Section 3, it is possible to have bounded functions with signed α-order Caputo derivative, but not converging as t → ∞, which is a surprising result interesting by itself.

To overcome this difficulty and obtain the convergence on Z, we need to give additional assumptions on the geometry of Z and argmin{g}. Precise assumptions are stated in Section 4 but basically, we assume

Z ∩ argmin{g} = ∅, (1.14)
for any z ∈ Z, there exists a rectifiable curve in Z joining argmin{g} and z. (1.15) These assumptions are inspired from the classical case, where it is known that the solution u of (1.8) is the value function of an optimal control problem for which assumptions (1.14)-(1.15) means roughly that we can travel with minimal cost on Z. But let us point out that these geometrical assumptions are not required in the classical case and that there is no rigorous link between the fractional case and the expected control problem, which could make these ideas rigorous. See Camilli, De Maio, Iacomini [START_REF] Camilli | A Hopf-Lax formula for Hamilton-Jacobi equations with Caputo time derivative[END_REF] for the precise statement of a related control problem which should be associated with the case α ∈ (0, 1), and some discussion in this direction.

It follows that we need to translate these ideas in the PDE framework building some suitable supersolutions, which tends to 0 thanks to (1.14)- (1.15). By comparison, we obtain estimates min

T N g ≤ u(z, t) + ct α ≤ min T N g + Lip(g) length(γ) E(t), for any z ∈ Z,
where γ is a rectifiable curve on Z joining argmin{g} and z, and E is the solution of a fractional ODE with limit 0 at +∞; see Theorem 4.2 and Theorem 4.4 for an extension to some possibly infinite length curves. This implies the convergence of u on Z from which we deduce easily (1.13), see Corollary 4.6. This approach is new but we think it does not provide optimal results. In particular, it relies too heavily on the geometry of Z, which is not the case in the classical approach described above for α = 1. To go further, one would need some quantitative estimates on how much nonpositive is ∂ α t (u(x, t) + ct α ), which seems difficult to obtain, even in the case α = 1.

The paper is organized as follows. We start by introducing precisely the Caputo fractional operator and recalling some useful properties. Then we establish some regularity estimates for the solution of (1.1) in Section 2. Section 3 is devoted to a counterexample showing that a bounded function with nonnegative Caputo derivative does not necessarily converges. Finally, some positive results for the large time behavior of the solution of (1.1) are proved in Section 4.

Notations and preliminaries. We will write the Caputo derivative using (1.4) (with cα = 1 for simplicity). More precisely, let φ : (0, +∞) → R. We extend φ to (-∞, 0) by setting φ(s) = φ(0) for s < 0 and define, when it exists, for every t > 0 and 0 < δ < t,

(1.16) ∂ α t φ(t) = t -∞ φ(t) -φ(s) |t -s| 1+α ds = ∂ α t [t -δ]φ(t) + ∂ α t [t -δ, t]φ(t),
where, for a < b ≤ t, we set

∂ α t [a]φ(t) := a -∞ φ(t) -φ(s) |t -s| 1+α ds, ∂ α t [a, b]φ(t) = b a φ(t) -φ(s) |t -s| 1+α ds.
Notice that ∂ α t φ(t) is well-defined as soon as φ ∈ L 1 loc (0, +∞) and φ is C 1 in a neighborhood of t. More about the functional formulation of this operator can be found in the references [START_REF] Kilbas | Theory and applications of fractional differential equations[END_REF][START_REF] Diethelm | The analysis of fractional differential equations[END_REF][START_REF] Gorenflo | Mittag-Leffler functions, related topics and applications[END_REF]. For the definition of viscosity solutions to (1.1), we refer to [START_REF] Barles | Second-order elliptic integro-differential equations: viscosity solutions' theory revisited[END_REF][START_REF] Topp | Existence and uniqueness for parabolic problems with Caputo time derivative[END_REF].

Notice that for all β > 0, there exists c α,β > 0 such that

∂ α t t β = c α,β t β-α
, for all t ≥ 0. [10, Appendix B] (1.17)

We introduce the Mittag-Leffler functions E α (z) of order α as in [START_REF] Diethelm | The analysis of fractional differential equations[END_REF]Definition 4.1]. Recall that E α is smooth on R, E α (0) = 1 and we have the following useful properties t ∈ [0, +∞) → E α (-t) is positive, convex and nonincreasing, [22, Section 6.2], (1.18)

∂ α t E α (λt α ) = λE α (λt α ), for all t > 0 and λ ∈ R, [10, Theorem 4.3], (1.19) 1 Γ(1 -α)t ≤ E α (-t) ≤ 1 Γ(1 + α) -1 t
for all t ≥ 0, [START_REF] Simon | Comparing Fréchet and positive stable laws[END_REF]Theorem 4]. (1.20) We will write x = x/|x| for every x ∈ R N \ {0}.

Regularity Estimates

We start with a regularity result in time for bounded solutions to (1.1). This result is a consequence of the Hölder estimates reported by Barles, Chasseigne and Imbert in [START_REF] Barles | Hölder continuity of solutions of second-order nonlinear elliptic integro-differential equations[END_REF].

Theorem 2.1. Assume (1.5) or (1.6), and g ∈ Lip(T N ). Let u be bounded, continuous viscosity solution to (1.1)-(1.2). Then, u is Hölder continuous in time, that is there exists a constant L > 1 large enough such that

(2.1) |u(x, s) -u(x, t)| ≤ L|s -t| α , s, t ≥ 0.
The constant L depends on the data and ||u|| ∞ but does not depend on t.

Proof. Despite the proof we present here is valid for both (1.5) and (1.6), we underline the arguments in the later case. Let µ ∈ (0, 1] be away from zero, and denote ū = µu. If (1.6) is assumed, then we consider µ < 1, and if (1.5) is assumed, then µ = 1 in what follows.

Using the linearity of the Caputo derivative, it is direct to check that ū solves (2.2)

∂ α t ū + µH(x, µ -1 Dū) = 0 in Q, with initial condition ū(•, 0) = µg.
By (1.17) and g ∈ Lip(T N ), we have that g(x)±Ct α are respectively super-and subsolutions of (1.1)-(1.2) for C > 0 large enough depending only on g and H. Therefore, a direct application of the comparison results in [START_REF] Topp | Existence and uniqueness for parabolic problems with Caputo time derivative[END_REF] leads to the estimates

(2.3) g(x) -Ct α ≤ u(x, t) ≤ g(x) + Ct α , for all t ≥ 0.
These bounds can be readily adapted to ū by multiplying by µ the last inequality. By contradiction, assume that for all L > 1 we have sup

x∈T N ,s,t≥0 {u(x, s) -u(x, t) -L|s -t| α } =: θ L > 0. (2.4) Taking µ very close to 1 in term of θ L above, we can get sup x∈T N ,s,t≥0 {ū(x, s) -u(x, t) -L|s -t| α } = θ L /2 > 0.
For localization purposes, we introduce a function ψ β with the following properties: we consider ψ : R → R smooth and nondecreasing, such that

ψ(t) = 0 if t ≤ 1, ψ(t) ≥ 2||u|| ∞ if t ≥ 2,
and for β > 0 small we denote ψ β (t) = ψ(βt). Then, for all β small enough, we have max

x∈R N ,s,t∈R {ū(x, s) -u(x, t) -L|s -t| α -ψ β (s)} = θ L /4 > 0,
where we recall that we extend u(x, t) as g for negative times t.

Then, for β, > 0 small we define

Φ(x, y, s, t) := ū(x, s) -u(y, t) -L|s -t| α --2 |x -y| 2 -ψ β (s).
We have that Φ attains its maximum at a point (x, ȳ, s, t) ∈ Q2 and this maximum is bigger than θ L /4.

Standard arguments lead to

|x -ȳ| ≤ ω β ( ), |s|, | t| ≤ 2/β,
where

ω β ( ) → 0 as → 0 if β is fixed (ω β is a modulus of continuity in space of u in the compact set T N × [0, 2/β]), and
|s -t| ≤ C 0 L -1/α , (2.5)
for some constant C 0 > 0 just depending on ||u|| ∞ . In particular, for L large enough we may assume that |s -t| < 1.

Below, we use a constant C, which may vary line to line but only depend on the data of the problem and not on , β nor µ.

Here we claim that s, t > 0 for all L large, and small in terms of L. In fact, if s = 0 (the case t = 0 is analogous), then, using (2.3), we see that

θ L /4 < Φ(x, ȳ, s, t) ≤ µg(x) -g(ȳ) + (C 0 -L) tα ≤ C 0 L 0 ω β ( ) + C 0 (1 -µ) + (C 0 -L) tα ,
where L 0 is the Lipschitz constant of g. Taking L ≥ C 0 , we arrive at

θ L /4 ≤ C + C(1 -µ),
which is a contradiction if is taken small, and µ close to 1 in terms of L.

In addition, since u is uniformly continuous in space, then s = t for all L large enough and small in terms of L. Indeed, if s = t, then we would have

0 < θ L /4 ≤ ū(x, t) -u(ȳ, t) ≤ ω β ( ) + C 0 (1 -µ).
Hence, for L and β fixed, taking small and µ close to 1, we arrive at a contradiction. Since s, t > 0 and s = t, we can use the penalization defining Φ as a testing for u. For this, in what follows we write φ(x, y, s, t) := L|s -t| α + -2 |x -y| 2 + ψ β (s), from which we see that Φ(x, y, s, t) = ū(x, s) -u(y, t) -φ(x, y, s, t).

At this point, denoting φ 1 (x, s) = φ(x, ȳ, s, t) we notice that the function

(x, s) → ū(x, s) -(u(ȳ, t) + φ 1 (x, s))
has a local maximum point at (x, s), from which we can use the subsolution's viscosity inequality to write, for all δ > 0, that

∂ α t [s -δ]ū(x, •)(s) + ∂ α t [s -δ, s]φ 1 (x, •)(s) + µH(x, µ -1 Dφ 1 (x, s)) ≤ 0. Similarly, denoting φ 2 (y, t) = φ(x, y, s, t) we notice the function (y, t) → u(y, t) -(ū(x, s) -φ 2 (y, t))
has a local minimum point at (ȳ, t), from which we can use the supersolution's viscosity inequality to write, for all δ > 0, that

∂ α t [ t -δ]u(ȳ, •)( t) + ∂ α t [ t -δ, t](-φ 2 )(ȳ, •)( t) + H(ȳ, D(-φ 2 )(ȳ, t)) ≥ 0. Then, we subtract both inequalities to arrive at D ≤ H, (2.6) where, noticing that D(-φ 2 )(ȳ, t)-= Dφ 1 (x, s) = 2(x -ȳ)/ 2 , we write D =∂ α t [s -δ]ū(x, •)(s) + ∂ α t [s -δ, s]φ 1 (x, •)(s) -∂ α t [ t -δ]u(ȳ, •)( t) + ∂ α t [ t -δ, t]φ 2 (ȳ, •)( t). H =H(ȳ, p) -µH(x, µ -1 p), p := 2(x -ȳ)/ 2 .
From (1.6), we get

H ≤ -(1 -µ)c -1 H |p| m + c H (1 + |p| m ) + A(1 -µ) and from here, taking ≤ (1 -µ)c -2 H , we conclude that H ≤ C + A(1 -µ), (2.7)
where C > 0 just depends on c H . Consider δ = |s -t|/2 > 0 and we split the term D as

D = D 1 + D 2 with D 1 =∂ α t [s -δ]ū(x, •)(s) -∂ α t [ t -δ]u(ȳ, •)( t) D 2 =∂ α t [s -δ, s]φ 1 (x, •)(s) + ∂ α t [ t -δ, t]φ 2 (ȳ, •)( t).
We start with D 2 . Directly from the definitions, we see that

D 2 =L s s-δ |s -t| α -|z -t| α |s -z| 1+α dz + L t t-δ |s -t| α -|s -z| α | t -z| 1+α dz + s s-δ ψ β (s) -ψ β (z) |s -z| 1+α dz =:D 2 + o β (1),
where o β (1) → 0 as β → 0 independent of the rest of the variables by the smoothness of ψ. Now, performing the change of variables z = s + y in the first integral, and z = t + y in the second, we arrive at

D 2 = L 0 -δ |s -t| α -|s -t + y| α |y| 1+α dy + L 0 -δ |s -t| α -|s -t -y| α |y| 1+α dy = L 0 -δ |s -t| α -|s -t + y| α |y| 1+α dy + L δ 0 |s -t| α -|s -t + y| α |y| 1+α dy = -L δ -δ |s -t + y| α -|s -t| α |y| 1+α dy = -L δ -δ |s -t + y| α -|s -t| α -α|s -t| α-1 s -t y |y| 1+α dy,
where the last equality comes from the symmetry of the kernel. Performing a second order expansion and recalling that δ = |s -t|/2, we obtain that there exists ρ(y) ∈ (-δ, δ) such that

-|s -t + y| α -|s -t| α -α|s -t| α-1 s -t y = - α(α -1) 2 |s -t + ρ(y)| α-2 y 2 ≥ - α(α -1) 2 α-1 |s -t| α-2 y 2 .
Therefore,

D 2 ≥ - Lα(α -1) 2 α-1 |s -t| α-2 δ -δ |y| 1-α dy = α(1 -α) 2 -α L =: cL,
and from here

D 2 ≥ cL -o β (1). (2.8)
Now we deal with D 1 , writing

D 1 =∂ α t [s -1]ū(x, •)(s) -∂ α t [ t -1]u(ȳ, •)( t) + ∂ α t [s -1, s -δ]ū(x, •)(s) -∂ α t [ t -1, t -δ]u(ȳ, •)( t) ≥ -C||u|| ∞ + ∂ α t [s -1, s -δ]ū(x, •)(s) -∂ α t [ t -1, t -δ]u(ȳ, •)( t) =: -C||u|| ∞ + D 1 ,
where C > 0 depends only on α.

To deal with the remaining term D 1 , we assume that s < t (the case t < s follows the same lines). Performing similar change of variables as above and using the maximal inequality Φ(x, ȳ, s, t) ≥ Φ(x, ȳ, s + y, t + y) we arrive at

D 1 = -δ -1 ū(x, s) -u(ȳ, t) -(ū(x, s + y) -u(ȳ, t + y)) |y| 1+α dy ≥ -δ -1 ψ β (s + y) -ψ β (s) |y| 1+α dy.
Noticing that the smooth function

ψ β satisfies |ψ β | ≤ Cβ, we conclude that D 1 ≥ -o β (1).
From this

D 1 ≥ -C||u|| ∞ -o β (1).
Joining this with (2.7) and (2.8) in (2.6) we get

cL ≤ C||u|| ∞ + o β (1) + C + A(1 -µ).
Then, we let → 0 first, then µ → 1 and finally β → 0 and, having taken L large enough just in terms of the data and ||u|| ∞ , we reach a contradiction. It ends the proof. Now we would like to obtain estimates in space.

Theorem 2.2. Assume hypotheses of Theorem 2.1 hold. For each bounded viscosity solution to (1.1), there exists a modulus of continuity m ∈ C([0, +∞)) independent of t such that

|u(x, t) -u(y, t)| ≤ m(|x -y|) for all x, y ∈ T N , t ≥ 0.
If, in addition, H satisfies (1.7) for some m ≥ 1, then, for each β ∈ (0, 1), there exists a constant C > 0 such that

|u(x, t) -u(y, t)| ≤ C|x -y| β for all x, y ∈ T N , t ≥ 0.
The modulus m and the constant C depend on the datas and ||u|| ∞ but do not depend on t.

Proof. For ∈ (0, 1), we introduce the sup-convolution

u (x, t) = sup s≥0 {u(x, s) --1 |s -t| 2 }.
We collect some properties of this regularization of u:

(i) u is still Hölder continuous in time satisfying (2.1) like u,

(ii) ||u -u|| ∞ ≤ C α/2 , (iii) u is Lipschitz continuous in time with Lipschitz constant C -1 with C > 0 just
depending on ||u|| ∞ , (iv) u is a viscosity subsolution to (1.1) in T N × (a , +∞) for some a → 0 as → 0.

The proofs of (i) and (ii) are easy consequences of Theorem 2.1, (iii) is a classical property of the sup-convolution regularization and (iv) is proved in [START_REF] Topp | Existence and uniqueness for parabolic problems with Caputo time derivative[END_REF]. Now we prove the desired regularity of u by adapting the standard viscosity procedure to get regularity estimates from the coercivity of H. For any x 0 ∈ T N , s 0 > 0, β > 0 and L > 0 to be chosen, we consider (2.9) sup

x∈T N ,s≥0 {u (x, s) -u (x 0 , s 0 ) -L|x -x 0 | - |s -s 0 | 2 β 2 },
where is chosen small enough in order that s 0 > a in (iv). Classical results imply that this maximum is achieved at (x, s) with s → s 0 as β → 0. We take β small enough in order that s > a .

If x = x 0 , then we use (x, s) → u (x 0 , s 0 ) + L|x -

x 0 | + |s-s 0 | 2 β 2
as a test function for the subsolution u of (1.1) at (x, s) to get that, for every δ ∈ (0, 1), (2.10)

∂ α t [s -δ]u (x, •)(s) + ∂ α t [s -δ, s] | • -s 0 | 2 β 2 (s) + H(x, L x -x 0 ) ≤ 0. Actually, since |•-s 0 | 2 β 2
is smooth and u is Lipschitz continuous, we can send δ → 0 in the previous inequality. In other words, due to the Lipschitz continuity of u , we can use u itself as a test-function in the fractional derivative in the viscosity inequality, see [START_REF] Topp | Existence and uniqueness for parabolic problems with Caputo time derivative[END_REF]Proposition 2.4] for details.

It follows that it is enough to estimate the fractional term ∂ α t u (x, •)(s) that we expand, for δ > 0, as

∂ α t u (x, •)(s) = ∂ α t [s -1]u (x, •)(s) + ∂ α t [s -1, s -δ]u (x, •)(s) + ∂ α t [s -δ, s]u (x, •)(s). At first, from (ii), ∂ α t [s -1]u (x, •)(s) ≥ -C||u || ∞ ≥ -C(||u|| ∞ + α/2
). Then, using (i),

∂ α t [s -1, s -δ]u (x, •)(s) = s-δ s-1 u (x, s) -u (x, s) |s -s| 1+α ds ≥ -C s-δ s-1 1 |s -s| ds ≥ -C| log(δ)|.
For the third term, we use (iii) to obtain

∂ α t [s -δ, s]u (x, •)(s) ≥ -C s s-δ 1 |s -s| α ds ≥ - C δ 1-α .
Plugging these estimates in (2.10), we obtain

H(x, L x -x 0 ) ≤ C(1 + | log(δ)| + -1 δ 1-α ,
where C > 0 just depends on the data and ||u|| ∞ . Taking the minimum for δ > 0 we arrive at

(2.11) H(x, L x -x 0 ) ≤ C(1 + | log( )|).
From the coercivity of H, we reach a contradiction if L = L( ) is large enough.

It follows that the maximum in (2.9) is achieved for x = x 0 , which implies

u (x, s) -u (x 0 , s) ≤ L( )|x -x 0 |, for all x ∈ T N .
Sending β → 0, recalling s → s 0 as β → 0 and that x 0 , s 0 are arbitrary, we finally obtain that, for all x, y ∈ T N , t > 0,

u(x, t) -u(y, t) ≤ u (x, t) -u (y, t) + C α/2 ≤ L( )|x -y| + C α/2 , 0 < < 1. (2.12)
This latter inequality means that u is uniformly continuous with respect to x independently of t.

In addition, if H satisfies (1.7) for some m ≥ 1, then (2.11) and (2.12) lead to

u(x, t) -u(y, t) ≤ C(1 + | log( )|)|x -y| + C α/2 , 0 < < 1.
Thus, taking the infimum with respect to 0 < < 1, we conclude that

|u(x, t) -u(y, t)| ≤ C(1 + | log |x -y||)|x -y|,
from which the result follows.

Remark 2.3. When H has a sublinear growth, it is possible to obtain some better regularity estimates in space, namely, Lipschitz estimates. More precisely, if (1.5) holds and g ∈ Lip(T N ), then every bounded viscosity solution to (1.1) satisfies

|u(x, t) -u(y, t)| ≤ (1 + Lip(g)) E α (2c H t α ) |x -y| for all x, y ∈ T N , t ≥ 0.
Such a result was already showed in Giga and Namba [START_REF] Giga | Well-posedness of Hamilton-Jacobi equations with Caputo's time fractional derivative[END_REF]. We do not focus on such results because the Lipschitz constant depends heavily on time, a dependence we want to avoid in order to obtain the large time behavior of the solution.

3. Oscillating function with positive Caputo derivative.

In this section we construct a bounded function u

: [0, +∞) → R such that ∂ α t u ≥ 0 but such that lim inf t→+∞ u(t) < lim sup t→+∞ u(t),
which prevents u to have any limit as t → +∞. This result shows a great contrast with the standard case α = 1 in which ∂ t u ≥ 0 implies that u is a nondecreasing function and therefore it is convergent.

In what follows, for any α ∈ (0, 1), we define the incomplete regularized beta function (see [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]Chapter 6]) by

B α [z 0 , z 1 ] := 1 π csc(απ) z 1 z 0 t -α (1 -t) α-1 dt, for all 0 ≤ z 0 ≤ z 1 ≤ 1,
and we simply write

B α [z] = B α [0, z]. We remark that B α [0, 1] = 1 ([1, 6.1.

and 6.2.2]).

We also define

b α := B -1 α [1/2] ∈ (0, 1), where B -1 α [•] is the inverse function of B α [•].
As an example, if α = 1/2, then b α = 1/2. Hence, in the general case, with this choice of b α , we have

B α [0, b α ] = B α [b α , 1] = 1/2. (3.1)
Then, we define

η α := π csc(απ)B α [b 3 α , b 2
α ] ∈ (0, 1), and consider the continuous functions

f 1 (t) = 1 if t ∈ [0, 1], t -α if t ≥ 1,
and

f 2 (t) =        t-a 2k k if t ∈ [a 2k , a 2k + k ), 1 if t ∈ [a 2k + k , a 2k+1 -k ), a 2k+1 -t k if t ∈ [a 2k+1 -k , a 2k+1 ), 0
if not, where

a k := (1/b α ) k and k := 1 -b 2 α 4 η α a 2k
, for all k ≥ 0.

Next, we consider the continuous function f := f 1 f 2 (see Figure 1) and we define u as

(3.2) u(t) = t 0 f (z) (t -z) 1-α dz, for all t ≥ 0. Figure 1. Behavior of f = f 1 f 2 t 1 t α f a 0 a 1 a 2 a 3 a 4 a 5
The function f is regular enough (locally Lipschitz) to use the representation formula in [10, Theorem 3.7], meaning that u solves the fractional ODE

(3.3) ∂ α t u = f in (0, +∞), u(0) = 0,
where ∂ α t is the Caputo derivative of order α ∈ (0, 1). Notice that u = 0 in [0, 1] and 0 ≤ u is bounded in R + . In fact, for t > 1 we see that

u(t) ≤ t 1 f 1 (z) (t -z) 1-α dz, from which we get that u(t) ≤ t α-1 t 1 z -α (1 -z/t) 1-α dz = t α-1 1 1/t t -α y -α (1 -y) 1-α tdy ≤ 1 0 dy y α (1 -y) 1-α = π csc(απ),
from which u is bounded. Now we compare u(a 2N +1 ) and u(a 2N +2 ), for N large enough such that (3.4)

1 1-N /a 2N +1 dy y α (1 -y) 1-α < η α /4 and a 2N (a 1 -1) = (1/b α ) 2N (1/b α -1) ≥ 2.
Using the definition of f we see that

u(a 2N +1 ) = a 2N +1 1 z -α f 2 (z) (a 2N +1 -z) 1-α dz = N k=0 a 2k + k a 2k z -α f 2 (z) (a 2N +1 -z) 1-α dz + a 2k+1 -k a 2k + k z -α (a 2N +1 -z) 1-α dz + a 2k+1 a 2k+1 -k z -α f 2 (z) (a 2N +1 -z) 1-α dz = N k=0 a 2k+1 a 2k z -α (a 2N +1 -z) 1-α dz - N k=0 a 2k + k a 2k z -α (1 -f 2 (z)) (a 2N +1 -z) 1-α dz - N k=0 a 2k+1 a 2k+1 -k z -α (1 -f 2 (z)) (a 2N +1 -z) 1-α dz =:v 1 (a 2N +1 ) -v 2 (a 2N +1 ) -v 3 (a 2N +1 ),
and similarly

u(a 2N +2 ) =: v 1 (a 2N +2 ) -v 2 (a 2N +2 ) -v 3 (a 2N +2 ).
From here, by simple integration we get

v 1 (a 2N +1 ) = N k=0 a 2k+1 /a 2N +1 a 2k /a 2N +1 dy y α (1 -y) 1-α =π csc(απ) N k=0 B α [a 2k /a 2N +1 , a 2k+1 /a 2N +1 ]. Moreover, v 1 (a 2N +2 ) = π csc(απ) N k=0 B α [a 2k /a 2N +2 , a 2k+1 /a 2N +2 ].
Now, we estimate the term v 1 (a 2N +2 ) -v 1 (a 2N +1 ). For this, using the definition of b α , η α and (3.1), we notice that

v 1 (a 2N +2 ) -v 1 (a 2N +1 ) π csc(απ) = N k=0 (B α [a 2k /a 2N +2 , a 2k+1 /a 2N +2 ] -B α [a 2k /a 2N +1 , a 2k+1 /a 2N +1 ]) ≤B α [0, b α ] - N k=0 B α [a 2k /a 2N +1 , a 2k+1 /a 2N +1 ] =B α [0, b α ] -B α [b α , 1] - N -1 k=0 B α [a 2k /a 2N +1 , a 2k+1 /a 2N +1 ] = - N -1 k=0 B α [a 2k /a 2N +1 , a 2k+1 /a 2N +1 ] ≤ -B α [b 3 α , b 2 α ]
. Therefore, we have that v 1 (a 2N +2 ) -v 1 (a 2N +1 ) ≤ -η α . Hence, using the above result and the fact that v 2 , v 3 ≥ 0, we have that

u(a 2N +2 ) -u(a 2N +1 ) ≤ -η α + v 2 (a 2N +1 ) + v 3 (a 2N +1 ).
Finally, we estimate the last two terms

v 2 (a 2N +1 ) = N k=0 a 2k + k a 2k z -α (1 -z-a 2k k ) (a 2N +1 -z) 1-α dz = N k=0 k 1 0 1 -y (a 2k + k y) α (a 2N +1 -a 2k -k y) 1-α dy.
We have a 2k + k y ≥ a 0 = 1 and, using k ≤ 1 and (3.4),

a 2N +1 -a 2k -k y ≥ a 2N +1 -a 2N -1 = a 2N (a 1 -1) -1 ≥ 1. It follows v 2 (a 2N +1 ) ≤ N k=0 k 1 0 dy = ∞ k=0 k = η α 4 .
Similarly, we have that

v 3 (a 2N +1 ) ≤ N -1 k=0 a 2k+1 a 2k+1 -k z -α (1 - a 2k+1 -z k ) (a 2N +1 -z) 1-α dz + a 2N +1 a 2N +1 -N z -α (a 2N +1 -z) 1-α dz = N -1 k=0 k 1 0 1 -y (a 2k+1 -k y) α (a 2N +1 -a 2k+1 + k y) 1-α dy + a 2N +1 a 2N +1 -N z -α a 1-α 2N +1 (1 -z/a 2N +1 ) 1-α dz.
To estimate the first term above we notice first that a (3.4). To estimate the second term, we notice

2k+1 -k y -1 ≥ a 1 -(1 -b 2 α )/4 = (b α -1)(b 2 α -b α -4)/(4b α ) ≥ 0 since 0 < b α < 1. Moreover, for k ≤ N -1, a 2N +1 -a 2k+1 + k y ≥ a 2N +1 -a 2N -1 ≥ a 2N (a 1 -1) ≥ 2 by
a 2N +1 a 2N +1 -N z -α a 1-α 2N +1 (1 -z/a 2N +1 ) 1-α dz = 1 1-N /a 2N +1 1 y α (1 -y) 1-α dz ≤ η α 4 ,
using again (3.4). It follows

v 3 (a 2N +1 ) ≤ N -1 k=0 k 1 0 dy + η α 4 = η α 2 .
Therefore, u(a 2N +2 ) -u(a 2N +1 ) ≤ -η α /4 < 0, and this means that

lim inf t→∞ u(t) -lim sup t→∞ u(t) ≤ -η α /4,
from which u does not have any limit at infinity.

Ergodic Large Time Behavior.

In this section we present some cases for which ergodic large time behavior (1.13) holds. The main assumption here follows the classical requirements of Namah and Roquejoffre [START_REF] Namah | Remarks on the long time behaviour of the solutions of Hamilton-Jacobi equations[END_REF], see Assumptions (1.9).

We have in mind the classical Eikonal case

H(x, p) = F (x, p) -f (x) = a(x)|p| -f (x), (4.1)
where a, f : T N → R are Lipschitz continuous, a(x) ≥ a > 0 and f (x) ≥ min T N f = 0. However, we are able to deal with Hamiltonians with superlinear growth on the gradient.

We notice that the convexity and the coercivity condition in (1.5) leads to a quantitative growth for the Hamiltonian, that is

H(x, p) ≥ C -1 |p| -C, for all x ∈ T N , p ∈ R N ,
for some constants C > 0. Since a similar condition is found when (1.6) holds, throughout this section we assume the existence of a constants C H > 1 and m ≥ 1 such that

(4.2) F (x, p) ≥ C -1 H |p| m -C H , for all x ∈ T N , p ∈ R N .
We also require some assumption on the behavior of F near p = 0 which is not reflected by (4.2). In order to be able to deal with more general Hamiltonian, e.g., smooth ones which are nonnegative and nondegenerate near p = 0, we introduce an additional assumption:

There exists ν, r > 0 and k ≥ 1 such that F (x, p) ≥ ν|p| k for all x ∈ T N , p ∈ B(0, r). (4.3) Thus, if (1.5) or (1.6) holds, then the later condition together with (1.9) lead to for every R > 0, there exists ν R > 0 such that F (x, p) ≥ ν R |p| k for all p ∈ B(0, R). (4.4) This is a sort of nondegeneracy condition in the sense that F is not too flat around p = 0.

By replacing f with f -min T N f in (1.10), we may assume without loss of generality that min

T N f = 0.
It follows that c = 0 in (1.12).

As we will see later in the proof of Lemma 4.3, the solutions of (1.1) for x ∈ Z are strongly related to the solutions of the ODE ∂ α t E(t) + A|E(t)| k = 0, for which we state a technical lemma.

Lemma 4.1. For every A > 0 and k ≥ 1, there exists a unique positive solution

E ∈ C([0, ∞)) ∩ C 1 ((0, ∞)) to ∂ α t E(t) + A|E(t)| k = 0, E(0) = 1, (4.5)
such that E(t) 0 as t → +∞. Moreover, there exists C = C(A, k, α) > 0 and, for all > 0, there exists

C = C (A, k, α) such that C t α/k ≤ E(t) ≤ C t α/k-for t large enough. (4.6)
Proof of Lemma 4.1. Existence and uniqueness of the positive decreasing solution E of (4.5) satisfying the lower bound in (4.6) is given by [START_REF] Feng | Continuous and discrete one dimensional autonomous fractional ODEs[END_REF]Theorem 5.10].

Concerning the (upper) estimates in (4.6), we start with the case k = 1. In this case the related ODE (4.5) reads ∂ α t E(t) + AE(t) = 0, for which we have the explicit solution t → E α (-At α ), see (1.19). In this case, the estimate (4.6) is then optimal and given by (1.20). Now we concentrate on the case k > 1. Below, c is a positive constant which may change line to line. Also, by

∂ α t (1 + t) -p we mean ∂ α t (1 + •) -p (t).
Claim: For each p > 0, there exists c > 0 just depending on α and p such that

-ct -α ≤ ∂ t (1 + t) -p ≤ 0 for all t > 1.
The upper bound is obvious. The lower bound can be obtained by a combination of [10, p.193] and [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]Chapter 15,7.3], but we present here an alternative proof for completeness.

Using the definition of Caputo derivative, for t > 1 we see that

∂ α t (1 + t) -p ≥ ∂ α t [t/2, t](1 + t) -p + ∂ α t [0, t/2](1 + t) -p + 0 -∞ -1 |t -z| 1+α dz ≥ ∂ α t [t/2, t](1 + t) -p + ∂ α t [0, t/2](1 + t) -p -ct -α ,
for some c > 0 just depending on α. Using the Mean Value Theorem, there exists a constant c depending on p such that

∂ α t [t/2, t](1 + t) -p ≥ -c(1 + t) -p-1 t t/2 |t -z| -α dz ≥ -ct -p-α ,
meanwhile, neglecting positive terms, we see that

∂ α t [t/2, t](1 + t) -p ≥ - t/2 0 (1 + z) -p |t -z| 1+α dz ≥ - t/2 0 dz |t -z| 1+α ≥ -ct -α ,
and joining the above inequalities we conclude the Claim. Let > 0 be small enough in order to have p := α/k -> 0. Take C > 0 large enough such that C(1 + t) -p ≥ C2 -p ≥ E in [0, 1]. By the Claim, it is possible to take C larger if it is necessary to get

∂ α t C(1 + t) -α/k+ + (C(1 + t) -α/k+ ) k ≥ 0 in [1, +∞).
Then, by comparison, we arrive at 0 ≤ E(t) ≤ C (1 + t) -α/k+ for all t ≥ 0. This concludes the proof.

In order to state our key result to obtain the large time behavior, we need some definitions. Given two points x 0 , x 1 ∈ R N , we denote [x 0 , x 1 ] the line segment joining x 0 and x 1 . For a set of points x 0 , x 1 , ..., x n , with n ∈ N, we denote

[x 0 , ..., x n ] = n i=1 [x i , x i-1 ],
that is, the polygonal curve joining the points x i , i = 0, ..., n. The length of a finite polygonal

line [x 0 , x 1 , • • • , x n ], x i ∈ T N , is given by ([x 0 , x 1 , • • • , x n ]) = n i=1 |x i -x i-1 |. A continuous curve γ : [0, 1] → T N is said to be rectifiable if (γ) := sup n ∈ N 0 := t 0 < t 1 < • • • < tn = 1 ([γ(t 0 ), γ(t 1 ), • • • , γ(t n )]) < +∞.
We call (γ) the length of γ.

As explained in the introduction, the proof of Namah-Roquejoffre Theorem [20, Theorem 1] relies on three steps, the limiting one being the second one, i.e., to prove that u(•, t) converges on Z = {f = 0} when t → +∞ (recall that we assume c = 0). We prove now this key result under the additional assumption

argmin{g} ∩ Z = ∅. (4.7)
The large time behavior is an easy consequence, see Corollary 4.6. Theorem 4.2. Assume (1.5) or (1.6), (1.9), (4.3) and (4.7). Assume further that for each z ∈ Z, there exists x 0 ∈ Z ∩ argmin{g} and a rectifiable curve γ : [0, 1] → T N such that γ(0) = x 0 , γ(1) = z and γ(t) ∈ Z for all t ∈ [0, 1]. Then, the unique solution u to (1.1)-(1.2) converges on Z, i.e., for every x ∈ Z, u(x, t) → min{g} as t → +∞.

Before giving the proof of the theorem, we state the following key lemma. where E(t) 0 as t → +∞ is a function which depends on H, ||f || ∞ , Lip(g), N and (γ).

Proof of Lemma 4.3. Without loss of generality we can assume min{g} = 0. From this and (1.9), we obtain that 0 is a subsolution of (1.1)-(1.2) in Q. Therefore, by comparison, 0 ≤ u in Q.

For the upper bound, the idea is to construct a function U such that U (z, t) 0 as t → +∞ and such that u(z, t) ≤ U (z, t) for all t > 0. This is performed by an inductive procedure, building a sequence of functions (U i ) 0≤i≤n-1 which are supersolutions for the equation solved by u but in the set +∞), with some control on the line {x i } × [0, +∞) in order to use comparison principles for the Cauchy-Dirichlet problem.

Q i := Q \ {x i } × [0,
We divide the proof in several steps.

Step 1. Definition of E(t). By Lemma 4.1, for every A > 0 and k ≥ 1, there exists a unique positive solution

E ∈ C([0, ∞)) ∩ C 1 ((0, ∞)) to (4.5) such that E(t)
0 as t → +∞. Notice that, since E is nonincreasing, ∂ α t E(t) ≤ 0. We now define A and other constants, the definition of which will be clear below. We set L := Lip(g) and

M := L + C 2 H + C H ||f || ∞ + LC H ( √ N + (γ)), (4.9)
where C H appears in (4.2) and √ N = diam(T N ). From (4.3) and (4.4), we may define ν L+M > 0 such that

F (x, p) ≥ ν L+M |p| k for |p| ≤ L + M. (4.10)
We then fix

A = ν L+M L k-1 √ N + (γ) (4.11)
in (4.5). Notice that we may assume without loss of generality that that A ≤ 1 by decreasing ν L+M > 0 if necessary.

Step 2. Definition of the function U i , 0 ≤ i ≤ n -1. We set

U 0 (x, t) = L|x -x 0 |E(t) + M d [x 0 ,x 1 ] (x) and (4.12) U i (x, t) = L i j=1 |x j -x j-1 |E(t) + L|x -x i |E(t) + M d [x i ,x i+1 ] (x), 1 ≤ i ≤ n -1,
where M is given by (4.9) and d [x i ,x i+1 ] denotes the (periodic) distance function to the segment [x i , x i+1 ], that is, for each x ∈ T N (cast as a point in [0, 1) N ), we write

d [x i ,x i+1 ] (x) = inf y∈[x i ,x i+1 ],k∈Z N |x + k -y|.
This is a 1-Lipschitz continuous. At the points where it is differentiable we have the gradient meets x -p i , where p i is the projection of x to the segment, from which |Dd [x i ,x i+1 ] (x)| = 1. In addition, for each point in the set of non differentiability of d [x i ,x i+1 ] which do not lie on the segment [x i , x i+1 ], there is not C 1 function touching the function from below.

Step 3. The initial supersolution U 0 . We prove actually that U 0 is a supersolution. At first, for t = 0 and all x ∈ T N , since L = Lip(g) and g(x 0 ) = 0, we have

U 0 (x, 0) ≥ L|x -x 0 | ≥ g(x) -g(x 0 ) = g(x).
If t > 0, x ∈ [x 0 , x 1 ] and U 0 is C 1 at (x, t), by the choice of the constant M 0 we use the coercivity properties of the Hamiltonian to write the following computation holds

∂ α t U 0 + F (x, DU 0 ) -f (x) = L|x -x 0 |∂ α t E(t) + F x, LE(t) x -x 0 + M 0 Dd [x 0 ,x 1 ] (x) -f (x) ≥ L √ N ∂ α t E(t) + 1 C H LE(t) x -x 0 + M 0 Dd [x 0 ,x 1 ] (x) -C H -||f || ∞ ≥ M C H - L C H -C H -||f || ∞ + L √ N ∂ α t E(t) + A|E(t)| k -L √ N A|E(t)| k ≥ M C H - L C H -C H -||f || ∞ -L √ N A,
where we used (4.2), 0 ≤ E(t) ≤ 1, ∂ α t E(t) ≤ 0, (4.5) and

√ N = diam T N .
By the choice of M in (4.9) and since A ≤ 1, the right hand side of the previous inequality is nonnegative. So U 0 is a supersolution outside [x 0 , x 1 ]. Now, let t > 0, x ∈ [x 0 , x 1 ] and consider any C 1 test-function ϕ touching U 0 from below at (x, t), i.e., U 0 (y, s) ≥ ϕ(y, s) and U 0 (x, t) = ϕ(x, t). Since U 0 is sufficiently smooth in time, the following computation holds in a classical way,

∂ α t ϕ(x, t) = t -∞ ϕ(x, t) -ϕ(x, s) |t -s| 1+α ds ≥ t -∞ U 0 (x, t) -U 0 (x, s) |t -s| 1+α ds = ∂ α t U 0 (x, t) = L|x -x 0 |∂ α t E(t).
Since the right hand side of the above inequality is zero at x = x 0 and F (x, p) ≥ 0, the supersolution property holds at x = x 0 . It remains to prove that U 0 is a supersolution for x ∈ [x 0 , x 1 ] \ {x 0 }. In this case, x -

h x -x 0 ∈ [x 0 , x 1 ] for h > 0 enough small and d [x 0 ,x 1 ] (x) = d [x 0 ,x 1 ] (x -h x -x 0 ) = 0. It follows ϕ(x -h x -x 0 , t) -ϕ(x, t) ≤ U 0 (x -h x -x 0 , t) -U 0 (x, t) = L|x -h x -x 0 -x 0 |E(t) -L|x -x 0 |E(t) = -hE(t).
Dividing by h and letting h 0 we obtain a lower-bound for |Dϕ(x, t)|,

|Dϕ(x, t)| ≥ Dϕ(x, t), x -x 0 ≥ LE(t). (4.13)
Noticing that U 0 is Lipschitz continuous with constant L+M in space, we have also an upperbound |Dϕ(x, t)| ≤ L + M . For x ∈ [x 0 , x 1 ] \ {x 0 }, recalling f (x) = 0, we use the behavior of F near the origin in (4.4) together with (4.13) to write

∂ α t ϕ(x, t) + F (x, Dϕ(x, t)) -f (x) ≥ ∂ α t U 0 (x, t) + ν L+M |Dϕ(x, t)| k = L|x -x 0 |∂ α t E(t) + ν L+M L k |E(t)| k ≥ L √ N ∂ α t E(t) + A|E(t)| k -L √ N A|E(t)| k + ν L+M L k |E(t)| k ≥ L ν L+M L k-1 -A √ N |E(t)| k
by (4.5). Recalling the choices of A in (4.11), the right hand side of the above inequality is nonnegative.

Since the viscosity inequality follows at once in the points where the distance function cannot be touched from below, the previous discussion leads to conclude that U 0 is a supersolution of (1.1)-(1.2) in Q. By comparison, we obtain u(x, t) ≤ U 0 (x, t) for all (x, t) ∈ Q, from which, in particular we get that u(x 1 , t) ≤ L|x 1 -x 0 |E(t) for all t.

Step 4. Proof by induction for U i , i ≥ 1. By induction, we will prove that U i in (4.12) satisfies

∂ α i U i + H(x, DU i ) ≥ 0 in Q i , u ≤ U i in ∂ p Q i , (4.14)
where ∂ p Q i is the parabolic boundary {x i } × (0, t) ∪ T N × {0}.

We first deal with the Cauchy-Dirichlet condition. By assumption, we have 0 ≤ u(x i , t) ≤ U i-1 (x i , t) for every t ≥ 0. When i = 1, it follows

u(x 1 , t) ≤ U 0 (x 1 , t) = L|x 1 -x 0 |E(t) ≤ √ N LE(t) ≤ U 1 (x 1 , t). When i ≥ 2, we have u(x i , t) ≤ U i-1 (x i , t) = L i-1 j=1 |x j -x j-1 |E(t) + L|x i -x i-1 |E(t) = U i (x i , t).
For t = 0 and all x ∈ T N , we have, using the triangle inequality,

U i (x, 0) = L i j=1 |x j -x j-1 |E(0) + L|x -x i |E(0) ≥ L|x -x 0 | ≥ g(x) -g(x 0 ) = g(x),
so the initial condition is satisfied. Now we deal with the PDE in (4.14). The proof follows the same lines as the one in Step 4, so we only sketch it.

If t > 0 and x ∈ [x i , x i+1 ], then U i is C 1 at (x, t) and we can do the same computation as in Step 4 to obtain

∂ t α U i + F (x, DU i ) -f (x) = L (γ)∂ α t E(t) + L|x -x i |∂ α t E(t) + F x, L x -x i E(t) + M Dd [x i ,x i+1 ](x) -f (x) ≥ L( (γ) + √ N )∂ α t E(t) + M -L C H -C H -||f || ∞ ≥ L( (γ) + √ N ) ∂ α t E(t) + A|E(t)| k + M C H - L C H -C H -||f || ∞ -L( (γ) + √ N )A|E(t)| k .
By the choice of M in (4.9), recalling that A ≤ 1, we obtain that the right hand side of the above inequality is nonnegative. Now, let t > 0, x ∈ [x i , x i+1 ] \ {x i } and consider any C 1 test-function ϕ touching U i from below at (x, t), i.e., U i (y, s) ≥ ϕ(y, s) and U i (x, t) = ϕ(x, t). As in Step 4, we check easily that

∂ α t ϕ(x, t) ≥ ∂ α t U i (x, t) and, since d [x i ,x i+1 ] (x -h x -x i ) = 0 for h > 0 small enough, that |Dϕ(x, t)| ≥ LE(t). Moreover, |Dϕ(x, t)| ≤ L + M . It follows ∂ α t ϕ(x, t) + F (x, Dϕ(x, t)) -f (x) ≥ L( (γ) + √ N ) ∂ α t E(t) + A|E(t)| k -L( (γ) + √ N )A|E(t)| k + ν L+M L k |E(t)| k ≥ L ν L+M L k-1 -( (γ) + √ N )A |E(t)| k recalling that f = 0 on [x i , x i+1 ]
. By the choice of A in (4.11), the right hand side of the above inequality is nonnegative. It ends the proof of the supersolution property for U i , and therefore the inductive process. Then, using Cauchy-Dirichlet comparison principles in [START_REF] Namba | On existence and uniqueness of viscosity solutions for second order fully nonlinear PDEs with Caputo time fractional derivatives[END_REF] (or standard adaptation of the comparison principles presented in [START_REF] Topp | Existence and uniqueness for parabolic problems with Caputo time derivative[END_REF] when the Dirichlet condition is satisfied pointwisely), we conclude that u ≤ U i in Q i , from which we get u(z, t) ≤ U i (z, t) for all t > 0. This concludes the proof.

Proof of Theorem 4.2. We may assume without loss of generality that min g = 0. It follows that 0 ≤ u in Q and we need to prove that u(x, t) → 0 on Z as t → +∞.

Let > 0 and f := (f -C f ) + , where C f := Lip(f ). Then f is a periodic function that satisfies the following properties,

0 ≤ f ≤ f, ||f -f || ∞ ≤ C f , Z = {f = 0} ⊂ {d Z (x) ≤ } ⊂ Z = {f = 0}.
To prove the last inclusion, let x ∈ T N such that d Z (x) ≤ . Then there exists

x Z ∈ Z such that |x -x Z | ≤ and 0 ≤ f (x) ≤ f (x Z ) + C f |x -x Z | ≤ C f . Thus f (x) = 0. Now we consider (1.10)-(1.
2) by replacing f with f . Notice that the assumptions of Theorem 4.2 and Lemma 4.3 still holds. In particular, there exists a unique solution u . Moreover, u ± ||f -f ||c -1 α,α t α , where c α,α appears in (1.17), are respectively a super-and a subsolution of (1.10)-(1.2) with f . By comparison, we get

||u -u || ≤ C f c α,α t α . (4.15)
Let z ∈ Z and γ : [0, 1] → Z be a rectifiable polygonal curve such that γ(0) = x 0 ∈ Z ∩ argmin{g} and γ(1) = z. By assumption, there exists a sequence of subdivision t k 0 := 0 < t k 1 < • • • < t k n k := 1, k ∈ N and a finite polygonal line γ k := [γ(t k 0 ), • • • , γ(t k n k )] which satisfies (γ k ) (γ) as k → ∞. In particular, we can prove that dist(γ k , γ) → 0 as k → ∞. It follows that there exists k ∈ N such that γ k ⊂ Z and γ k is a finite polygonal line joining x 0 and z. We can apply Lemma 4.3 to obtain 0 ≤ u (z, t) ≤ Lip(g) (γ k )E(t), for all t ≥ 0. (4.16)

A priori, E depends on γ k through the dependence of A with respect to (γ k ) in (4.5) but, since (γ k ) ≤ (γ) < +∞, we can fix A and E independently of . From (4.15), we finally obtain 0 ≤ u(z, t) ≤ C f c α,α t α + Lip(g) (γ)E(t), for all t ≥ 0.

Sending first → 0 and then t → +∞, we conclude.

In the Eikonal case (4.1), we have an explicit formula for E in Lemma 4.3, see the proof of Lemma 4.1, which allows to deal with more involved Z. More precisely, we consider subsets Z satisfying the following assumption There exists D ≥ 1 and C > 0 such that, for all > 0 and x ∈ Z( ) := {d Z < }, there exists x 0 ∈ Z ∩ argmin{g} and a finite polygonal line γ ⊂ Z( ) such that γ is formed by at most C -D lines of length .

(4.17)

We remark that a set Z satisfying Assumption (4.17) is a curve of box-counting dimension D, see Falconer [START_REF] Falconer | Fractal geometry[END_REF]. In several interesting cases, box-counting dimension agrees with Hausdorff dimension and when D > 1 then the curve have infinite length. Proof of Theorem 4.4. We proceed exactly as in the proof of Theorem 4.2, where γ k is replacing with γ given by Assumption (4.17). From (4.15) and (4.16), we arrive at 0 ≤ u(z, t) ≤ C f c α,α t α + Lip(g) (γ )E(t), for all t ≥ 0. (4.18)

The difference with the proof of Theorem 4.2 is that γ is not necessarily rectifiable anymore. But we can estimate the length of γ thanks to (4.17) and take profit of the explicit formula for the solution of (4.5) in the Eikonal case k = 1.

The constant C > 0 below may change line to line but it does not depend neither on nor on t. We have (γ ) ≤ C -D , from Assumption (4.17 Plugging all the estimates in (4.18), we end up with 0 ≤ u(z, t) ≤ C t α + 1 2(D-1) t α , for all t ≥ 0, > 0. (4. [START_REF] Lions | Homogenization of Hamilton-Jacobi equations[END_REF] Minimizing over > 0, we obtain 0 ≤ u(z, t) ≤ Ct α 2D-3 2D-1 , and the right hand side tends to 0 as t → +∞ when D < 3 2 . Remark 4.5. Notice that the condition on D does not depend on α ∈ (0, 1). If we use the same approach in the classical case (α = 1), then we can repeat the above proof with E(t) = e -At and the right hand side of (4.19) now reads C t + 1-D e -t D-1 /C . We can prove that the minimum over > 0 tends to 0 as t → +∞ if and only if D < 2. Even, if we obtain a more general result than in the fractional case α ∈ (0, 1), this result is not optimal since the convergence on Z holds for any Z and argmin{g} even without any connectedness requirement (see Introduction).

Corollary 4.6. Under the assumptions of Theorem 4.2 or 4.4, the unique solution u of (1.1)-(1.2) satisfies u(x, t) + ct α -v(x) → 0 uniformly as t → +∞, where c = -min T N f and v is the unique solution of (1.12) satisfying v = min T N g on Z.

The proof of the corollary follows the procedure described in the introduction. We only sktech it. By comparison, we obtain that u + ct α is uniformly bounded in T N × [0, +∞). We then apply Theorems 2.1 and 2.2 to u + ct α to prove Step 1. Step 2 follows from Theorem 4.2 or 4.4 and Step 3 is classical. For details, we refer the reader to the survey of Barles in [START_REF] Achdou | Hamilton-Jacobi equations: approximations, numerical analysis and applications[END_REF] or Namah-Roquejoffre [START_REF] Namah | Remarks on the long time behaviour of the solutions of Hamilton-Jacobi equations[END_REF].

(1. 5 )

 5 |H(x, p) -H(y, p)| ≤ c H (1 + |p|)|x -y|, lim |p|→+∞ inf x∈T N H(x, p) = +∞.

  , µ -1 p) -H(x, p) ≥ (1 -µ) c -1 H |p| m -A , |H(x, p) -H(y, p + q)| ≤ c H (1 + |p| m )|x -y| + c H |q|(1 + |p| m-1 ),

Lemma 4 . 3 .

 43 Assume hypotheses of Theorem 4.3 hold. Let z ∈ Z, x 0 ∈ Z ∩ argmin{g} and assume that there exists a finite polygonal line γ := [x 0 , x 1 , • • • , x n := z] lying in Z and joining x 0 to z. Then, the unique solution u to (1.1)-(1.2) satisfies min{g} ≤ u(z, t) ≤ min{g} + Lip(g) (γ)E(t),(4.8)

Theorem 4 . 4 .

 44 (Eikonal case) Assume (1.5), (1.9), (4.4) with k = 1, and (4.17) for 1 ≤ D <

3 2 .

 2 Then, the unique solution u to (1.1)-(1.2) converges on Z, i.e., for every x ∈ Z, u(x, t) → min {g} as t → +∞.

1 C

 1 ),E(t) = E α (-At α ) ≤ C At α , from (1.19)-(1.20),A ≥ (γ ) , from (4.11) and Assumption (4.4) with k = 1.

Acknowledgements. Part of this work was done during a visit of E.T. and M.Y. to the Institut de Recherche Mathématique de Rennes. They acknowledge the hospitality of the Institut. O.L. is partially supported by the Agence Nationale de la Recherche (MFG project ANR-16-CE40-0015-01 and Centre Henri Lebesgue ANR-11-LABX-0020-01). E.T. was partially supported by Conicyt PIA Grant No. 79150056, Foncecyt Iniciación No. 11160817, and Dicyt -Apoyo Asistencia a Eventos 2018. M.Y. is partially supported by Escuela Politécnica Nacional, Proyecto PII-DM-2019-01.


Bertrand Plainchont

Patrick Giraudeau
email: patrick.giraudeau@univ-nantes.fr

Jean-Nicolas Dumez
email: jean-nicolas.dumez@univ-nantes.fr

Interleaved spatial/spectral encoding in ultrafast 2D NMR spectroscopy

Keywords: Ultrafast 2D NMR, spatial encoding, spatial/spectral pulses, excitation sidebands, interleaved excitation, TOCSY, HSQC

The possibility to record a full 2D spectrum in less than a second using ultrafast 2D NMR (UF2DNMR) is beneficial in many applications. However, the spatial encoding process on which UF2DNMR is based sets specific constraints on the spectral width and resolution of the resulting spectra. To overcome these limitations, a tailored encoding method using spatial/spectral pulses (SPSP) can be employed as an alternative to the traditional linear spatial encoding of interactions. Here we analyze and further develop this alternative spatial encoding strategy. We first carry out numerical simulations to describe the features of bidimensional SPSP pulses. Sidebands are identified along the spectral dimension of the excitation profile. An interleaved excitation scheme is then developed and implemented experimentally to suppress the unwanted signals that arise from these harmonic sidebands. Two examples are shown to illustrate the potential of the proposed approach. An ultrafast selective TOCSY spectrum is recorded to access sub-spectra and fully assign 1 H NMR resonances of individual residues of cyclosporin A. An ultrafast HSQC spectrum of a mixture of metabolites is recorded with an optimized spectral width in the spatially encoded dimension.

Introduction

Two-dimensional (2D) NMR spectroscopy is a rich source of information for structure and dynamics studies of a large variety of systems, with applications in chemistry, biology and medicine. The original and most frequently used scheme for the collection of 2D NMR data is based on an extension of Fourier transform NMR to a second dimension. The additional, indirect dimension is classically sampled in a discrete, stepwise manner, through the consecutive acquisition of 1D spectra where an additional amplitude or phase modulation of the signal is created as a function of an incremented time variable. 2D Fourier transformation of the acquired 2D time-domain data then provides a 2D frequency-domain dataset that can take multiple forms according to the large panoply of 2D correlation experiments [START_REF] Ernst | Principles of Nuclear Magnetic Resonance in One and Two Dimensions[END_REF]. However, this stepwise implementation results in long experiment durations, from several minutes to hours, which prevent the real-time observation of phenomena occurring in a short time scale. Many alternative methods have been developed with the motivation to perform 2D NMR experiments in a faster manner and unlock new applications or simply collect more 2D spectra in the same amount of time [START_REF] Rouger | Fast n-Dimensional Data Acquisition Methods[END_REF]. Ultrafast 2D NMR (UF2DNMR) is one of these methods and has been exploited during the last 15 years for a variety of applications such as reaction monitoring, complex mixture analysis or high-throughput metabolomics [START_REF] Gouilleux | Chapter Two -Ultrafast 2D NMR: Methods and Applications[END_REF].

UF2DNMR makes it possible to record a 2D NMR spectrum in a single scan of less than one second [START_REF] Frydman | The acquisition of multidimensional NMR spectra within a single scan[END_REF][START_REF] Frydman | Principles and Features of Single-Scan Two-Dimensional NMR Spectroscopy[END_REF]. In this approach, the indirect dimension is created by spatially encoding the spin interactions along one of the sample's dimension (typically z), i.e., imparting a phase variation of the form 𝜙(𝑧) = 𝐶 𝛺 ) 𝑧, where C is a spatial-encoding constant, to the spins' magnetization. A gradient-driven readout is used to decode the information. The resulting signal consists of a series of echoes, each corresponding to a different chemical-shift offset. In the most frequent implementation of UF2DNMR, spatial encoding is achieved with a pair of adiabatic frequency-swept pulses applied with bipolar gradient pulses, resulting in a linear relationship between the positions of the echoes in k-space (where k is the gradient area) and the chemical-shift values [START_REF] Pelupessy | Adiabatic Single Scan Two-Dimensional NMR Spectrocopy[END_REF]. The direct dimension is obtained with repeated bipolar gradient readout, analogous to echo planar spectroscopic imaging (EPSI) [START_REF] Mansfield | Multi-planar image formation using NMR spin echoes[END_REF], that allows the development of a signal modulation that is revealed by Fourier transform. The twodimensional data matrix acquired in the (k, t) domain yields a 2D spectrum in the (k, W) domain.

The speed of UF2DNMR experiments comes with a price. One limitation notably arises because of the use of pulsed-field gradients during acquisition: the large frequency dispersion induced by the magnetic field gradient requires a large receiver bandwidth, which results in a higher noise level than with conventional acquisition. In combination to possible diffusion effects, it causes important sensitivity penalty in the resulting spectra. This limitation can be counterbalanced with the coupling of UF2DNMR with hyperpolarization methods like dynamic nuclear polarization (DNP) [START_REF] Frydman | Ultrafast two-dimensional nuclear magnetic resonance spectroscopy of hyperpolarized solutions[END_REF][START_REF] Mishkovsky | Progress in Hyperpolarized Ultrafast 2D NMR Spectroscopy[END_REF] or parahydrogen-based signal amplification by reversible exchange (SABRE) [START_REF] Lloyd | Utilization of SABRE-Derived Hyperpolarization To Detect Low-Concentration Analytes via 1D and 2D NMR Methods[END_REF][START_REF] Daniele | Single-Scan Multidimensional NMR Analysis of Mixtures at Sub-Millimolar Concentrations by using SABRE Hyperpolarization[END_REF]. This is particularly relevant for single-shot polarization methods, which are incompatible with most classic 2D NMR pulse sequences. Another limitation is due to the capabilities of gradient hardware in high-resolution NMR, which limits the accessible spectral width and resolution in both dimensions for single-scan acquisitions. Several methodological developments that address the spectral width and resolution restrictions of UF2DNMR experiments have been reported [START_REF] Gouilleux | Chapter Two -Ultrafast 2D NMR: Methods and Applications[END_REF]. These include an interleaved acquisition method, in which multiple consecutive scans with shifted trajectories in the (k, t) domain are combined to increase the spectral width [START_REF] Frydman | Principles and Features of Single-Scan Two-Dimensional NMR Spectroscopy[END_REF]. A method to shift spectral regions in the spatially encoded dimension with a band-selective refocusing pulse flanked by folding gradients has been proposed by Pelupessy and coworkers [START_REF] Pelupessy | Improving resolution in single-scan 2D spectroscopy[END_REF]. An appropriate phase cycling of the selective pulse in multiple scans is used to unfold spectral regions. A so-called spatial/spectral (SPSP) encoding scheme has also been proposed by Shrot and Frydman [13]. This approach departs from the linear relation between frequency offset and echo position in k-space, and instead assigns arbitrary k-space positions to selected frequency offsets. It relies on the use of 2D selective pulses with both spatially and spectrally selective dimensions [START_REF] Meyer | Simultaneous spatial and spectral selective excitation[END_REF][START_REF] Bernstein | CHAPTER 5 -SPATIAL RADIOFREQUENCY PULSES[END_REF], initially developed for MRI applications to perform water/fat signal selection/elimination [START_REF] Block | Consistent fat suppression with compensated spectralspatial pulses[END_REF][START_REF] Schick | Simultaneous highly selective MR water and fat imaging using a simple new type of spectral-spatial excitation[END_REF][START_REF] Zur | Design of improved spectral-spatial pulses for routine clinical use[END_REF], and nowadays used extensively for hyperpolarized 13 C in vivo imaging [START_REF] Larson | Multiband excitation pulses for hyperpolarized 13C dynamic chemicalshift imaging[END_REF][START_REF] Larson | Investigation of tumor hyperpolarized [1-13C]-pyruvate dynamics using time-resolved multiband RF excitation echo-planar MRSI[END_REF][START_REF] Lau | Spectral-spatial excitation for rapid imaging of DNP compounds[END_REF]. While SPSP encoding alleviates some of the constraints that arise from the use of magnetic-field gradients during acquisition, it also comes with its own limitations. In particular, when a zig-zag trajectory is used in the (k, t) domain, the selection profile of SPSP pulses consists, in the spectral dimension, of a central band at the targeted offset, flanked by a series of harmonic excitation sidebands, as in the case of DANTE [START_REF] Morris | Selective excitation in Fourier transform nuclear magnetic resonance[END_REF] or binomial [START_REF] Hore | A new method for water suppression in the proton NMR spectra of aqueous solutions[END_REF][START_REF] Turner | Binomial solvent suppression[END_REF] pulses. A possible solution to avoid excitation of non-targeted resonances is to reject the sideband response towards empty spectral regions. Otherwise, without any control of the sidebands, the spectral information may be compromised by the presence of additional signals and may lead to wrong assignments and/or peak overlaps in 2D correlation spectra.

In the following sections, on the basis of a theoretical analysis and numerical simulations, we propose a method to eliminate undesirable excitation sidebands by interleaving trajectories in the excitation (k, t) domain. This approach is inspired by interleaved selection schemes developed in MRI to collect images or spectra of regions with a complex shape [START_REF] Hardy | 31P Spectroscopic localization using pinwheel NMR excitation pulses[END_REF][START_REF] Panych | Selection of high-definition 2D virtual profiles with multiple RF pulse excitations along interleaved echo-planar k-space trajectories[END_REF]. Applied to UF2DNMR, this method is well adapted to the study of samples with high dispersion of chemical shifts. A detailed description of SPSP pulse design with a focus on important parameters is provided to introduce interleaved SPSP excitation. The principle is then exemplified experimentally with both homonuclear and heteronuclear applications using two of the most commonly employed pulse sequences. A selective UF 2D TOCSY spectrum to access individual spin systems of cyclosporin A and an UF HSQC spectrum of a mixture of metabolites are recorded as illustrations.

Theory

Spatial encoding

UF2DNMR experiments include an excitation and encoding step followed by a mixing sequence and acquisition. Fig. 1a shows a pulse sequence scheme for the most frequent implementation of UF2DNMR, in which the encoding step consists of a pair of adiabatic frequency-swept pulses applied with bipolar gradient pulses. The effect of this block is the creation of a position-and frequency-offset-dependent phase for the transverse spin magnetization. This spatial encoding can be viewed as a spatial winding of the spin magnetization along the z-axis under a helical form (see Fig. 1b). Each individual helix has a specific pitch proportional to the chemical shift value. The unwinding of the spin magnetization during the gradient-driven acquisition leads to the detection of an echo for each encoded resonance frequency (Fig. 1c). The acquisition is most conveniently described by using the variable

𝑘(𝑡) = 𝛾 -. 𝐺 -(𝑡′)𝑑𝑡′ 2 3 (1)
that is commonly used in MRI, defining the so-called k-space. Echoes appear at positions ki in the k-space, the values of which have a linear relation with the frequency offsets Wi: 𝑘 i = -𝐶 𝛺 i , where C is a constant depending on the linear encoding method [START_REF] Tal | Single-scan multidimensional magnetic resonance[END_REF]. Assuming square shaped gradients, the span in k-space detectable during acquisition ∆𝑘 -depends on the acquisition gradient strength Ga and duration Ta:

∆𝑘 -= 𝛾 -𝐺 -𝑇 -. (2)
This k-space detection range defines the spectral width in the encoded dimension of the UF2DNMR spectra. Because of the linear relation between k and W, this spectral width SWk can be expressed in frequency units:

𝑆𝑊 : = ∆; < = = > < ? < @ < = . (3)
In this spatially encoded dimension, the detected echoes have a sinc shape, with a width of the order of L -1 , where L is the length of the sample detection area. In k-space, the resolution is proportional to the width of echoes. It comes that the resolution dnk, expressed in frequency units, is given by:

δ𝜈 : = C D= . (4)
In all linear encoding methods, the proportionality constant C depends on the encoding duration Te. For example, for the encoding scheme shown in Fig. 1a, 𝐶 = 𝑇 E /𝐿. Resolution can thus be improved by increasing the encoding duration but at the cost of sensitivity losses due to transverse relaxation and diffusion [START_REF] Giraudeau | Sources of sensitivity losses in ultrafast 2D NMR[END_REF][START_REF] Giraudeau | Resolution and sensitivity aspects of ultrafast J-resolved 2D NMR spectra[END_REF]. In the frequency dimension of the (k, W) domain, the spectral width SWW is simply given by the inverse of the sampling interval Δ𝑡 = 2𝑇 a :

𝑆𝑊 K = C L@ < . ( 5
)
After Fourier transform, any signal whose frequency lies outside the effective spectral width appears folded in the frequency domain in a similar way to what occurs in the indirect domain of conventional 2D NMR spectra when the Nyquist condition is violated. In UF2DNMR experiments with linear encoding, the spectral characteristics defined above are all closely related. A combination of Eqs. (3), ( 4) and ( 5) yields a relation giving an overview of the interdependence of the important parameters:

𝛾 -𝐺 -𝐿 = 2 MN O MN P QR O . ( 6
)
The relation between the two spectral widths in the (k, W) domain and the resolution in the spatially encoded dimension highlights the compromise that has to be made, which depends on the maximum available gradient strength. This is a key point for the setup of UF2DNMR experiments. In terms of spectral width, typical values for a single-scan experiment at 9.4 T on a conventional probe are 6 × 6 ppm for a 1 H-1 H COSY and 40 × 4 ppm for a 13 C-1 H HSQC [START_REF] Gouilleux | Ultrafast Multi-dimensional NMR: Principles and Recent Applications[END_REF]. The higher the magnetic field, the more problematic this limitation becomes. For example, at 16.4 T, single-scan spectral width values are reduced to 2.5 × 2.5 ppm for a 1 H-1 H COSY and 15 × 3 ppm for a 13 C-1 H HSQC.

Spatial/spectral encoding: concept

Spatial/spectral encoding (Fig. 1d) allows one to break away from the basic linear dependence between ki and Wi of the classic encoding scheme [START_REF] Shrot | Spatial/spectral encoding of the spin interactions in ultrafast multidimensional NMR[END_REF]. In this approach, multiple spin resonances are simultaneously manipulated to assign them the desired spatial windings (see Fig. 1e). In other words, the k-space position of the echo is arbitrarily chosen for each offset. This is achieved by using custom 2D spatial/spectral (SPSP) pulses to perform the encoding step (Fig. 1f). In this non-linear encoding scheme, the pseudo-frequency axis of the spatially encoded dimension does not exist anymore and signals are marked instead by k-axis indexes. The spectral width cannot be expressed in frequency units and the only relevant parameter is the maximum k-space range Dk. SPSP pulses permit to optimize the filling of the k-space with an evenly spacing of signals by a value dk at different indexes.

SPSP pulses involve the simultaneous application of radiofrequency pulses and oscillating gradients. As is usually the case for selective pulses, their design requires a priori knowledge of the targeted 1D resonance frequencies, which can be obtained from the analysis of a 1D spectrum. For UF2DNMR, the response of the transverse spin magnetization at the end of the encoding step corresponds to a target winding pattern for each selected resonance frequency with a specific pitch of the helix (equivalent to a k-space position). Once defined from a 1D spectrum analysis, (ki, Wi) pairs are used to build a target spatial/spectral profile for transverse spin magnetization with spatial winding in the (z, W) domain which is the sum of the individual desired response to the SPSP pulse:

𝑀 TU V-WXEY (𝑧, Ω) = \ 𝑒 ^;_ ` 𝑆(𝑧, Ω -Ω ) ) a )bC (7)
where 𝑆(𝑧, Ω -Ω ) ) is the local profile whose width in the frequency dimension defines the spectral resolution. Fig. 2a andb represent examples of the amplitude and phase components of a 2D target profile 𝑀 TU V-WXEY for three (ki, Wi) pairs. Here, the amplitude of the profile is set with a rectangular shape in spatial dimension and a Gaussian shape in spectral dimension for each of the three selected offsets (Fig 2a). The design of an SPSP pulse consists of finding the radiofrequency and gradient waveforms that yield, starting from equilibrium magnetization MZ, a transverse spin magnetization as close as possible to the target profile 𝑀 TU V-WXEY . Many approaches to the design of SPSP pulses have been proposed in the field of MRI [START_REF] Meyer | Simultaneous spatial and spectral selective excitation[END_REF][START_REF] Schick | Simultaneous highly selective MR water and fat imaging using a simple new type of spectral-spatial excitation[END_REF][START_REF] Zur | Design of improved spectral-spatial pulses for routine clinical use[END_REF][START_REF] Rieseberg | Two-dimensional spatially-selective RF excitation pulses in echo-planar imaging[END_REF][START_REF] Grissom | Maximum linear-phase spectral-spatial radiofrequency pulses for fat-suppressed proton resonance frequency-shift MR Thermometry[END_REF][START_REF] Schulte | Direct design of 2D RF pulses using matrix inversion[END_REF]. Here we follow the approach of Refs [START_REF] Shrot | Spatial/spectral encoding of the spin interactions in ultrafast multidimensional NMR[END_REF][START_REF] Schulte | Direct design of 2D RF pulses using matrix inversion[END_REF], by choosing a priori an encoding gradient waveform Ge(t), and relying on the small-tip-angle approximation [START_REF] Pauly | A k-space analysis of small-tip-angle excitation[END_REF] to obtain B1(t) from 𝑀 TU V-WXEY . In this approximation, the transverse spin magnetization created during a total pulse duration Tpulse is:

𝑀 TU (𝑧, Ω) ∝ . 𝐵 C (𝑡) 𝑒 ^K(@ efghi j2) 𝑒 ^;(2)@ efghi 3 𝑑𝑡 (8)
where

𝑘(𝑡) = ∫ 𝛾 E 𝐺 E (𝑡 l ) 𝑑𝑡′ @ efghi 2
is the k-space trajectory (Fig. 2c) defined by the integral of the encoding gradient shape Ge (Fig. 2e). Eq. ( 8) describes a 2D Fourier relationship between the transverse magnetization at the end of the pulse and the RF waveform when seen in the (k, t) domain. In analogy with the description of EPSI as a walk in detection (k, t) domain, SPSP pulses can be seen as a walk in excitation (k, t) domain.

The B1(t) radiofrequency shape for SPSP pulses is calculated by solving an expression where the target pattern ( 7) is substituted into Eq. ( 8). In practice, the equation is discretized on a grid that corresponds to the (k, t) trajectory for the selected gradient waveform. Since the grid is not Cartesian, a non-uniform Fourier transform (nuFT) algorithm is used to retrieve B1(t). In the example presented in Fig. 2, the radiofrequency shape (Fig. 2d) was calculated from the target profile with a nuFT algorithm based on conjugate gradients which converged within 3 iterations.

Spatial/spectral encoding: parameter relations

The relations between encoding/decoding parameters and the spectra's properties are modified when SPSP encoding is used, and depend on the choice of SPSP pulses. Here we choose to use a gradient waveform that consists of a train of bipolar square shaped pulses, shown in Fig. 2e, resulting in a zig-zag trajectory in excitation (k, t) domain, as shown in Fig. 2c. The 2D selective pulse then consists of a train of spatially selective subpulses, corresponding to each gradient lobe, and spectral selectivity is imposed by the global envelope of the pulse.

The total area of each encoding gradient lobe defines the span of k-space positions accessible during the encoding step. This is an important feature because signals acquired during detection must be resolved in the k-space dimension and thus separated by a minimum value dk (of the order of L -1 ). For m resonance frequencies to encode, the necessary k-space range is: Δ𝑘 E = 𝑚 δ𝑘. Thus, the encoding gradient strength Ge and oscillation period t should be adjusted according to the relation: Δ𝑘 E < 𝛾 E 𝐺 E 𝜏/2, assuming square shaped gradient lobes. Here, contrary to the linear encoding, the major limiting factor is the encoding gradient strength available to encode a given number of resonance frequencies with a given resolution.

As for linear encoding, for SPSP encoding the k-space range accessible during detection ∆𝑘 -and the spectral width in the frequency dimension SWW are given by Eqs. ( 2) and ( 5), respectively. Compared to linear encoding, if the filling of the k-space can be optimized with adequate resolution for the same number of signals, the acquisition gradient requires less strength. Alternatively, the same acquisition gradient can be used to detect a larger number of resolved signals in a given k-space range. Indeed, the role of SPSP pulses is akin to a controlled folding of signals in empty regions of the detectable k-space range. During acquisition, m signals with a width of the order of L -1 can be detected given the relation: ∆𝑘 -= 𝑚 𝐿 jC . Fig. 2. Spatial/spectral pulse design example. Three frequency offsets and the associated k-space positions define the target transverse spin magnetization profile 𝑀 TU V-WXEY in the (z, W) domain. The amplitude profile (a) is set with Gaussian shape in the spectral dimension and rectangular shape in the spatial dimension and the phase profile (b) is set according to Eq. ( 7). The k-space trajectory (c) corresponds to the integral of the encoding gradient Ge (e) whose oscillation period t is 1 ms and strength is 0.06 T.m -1 . A number of 25 zig-zag oscillating cycles are used to build the gradient shape leading to a total SPSP pulse duration of 25 ms. The radiofrequency shape (d) is calculated by nonuniform Fourier transform from the (z, W) to the (k, t) domain.

The encoding gradient oscillation period t plays a central role in SPSP encoding. The type of pulse chosen here, which consists of a train of small-tip-angle subpulses separated by delays, shares properties of DANTE [START_REF] Morris | Selective excitation in Fourier transform nuclear magnetic resonance[END_REF] or binomial [START_REF] Hore | A new method for water suppression in the proton NMR spectra of aqueous solutions[END_REF][START_REF] Turner | Binomial solvent suppression[END_REF] sequences. Specifically, the use of a discrete train of small-angle rotation, instead of a continuous rotation, results in the excitation of harmonic sidebands in addition to the targeted spectral region. Numerical spin simulations are helpful to characterize the excitation profile of SPSP pulses and in this case to detect and understand the occurrence of excitation sidebands. Fig. 3f shows the simulated transverse spin magnetization response |MXY(z, W)| to a single SPSP excitation (Fig. 3b), associated to a given k-space trajectory (Fig. 3d), and designed from the target profile 𝑀 TU V-WXEY depicted in Fig. 3a to select on-resonance frequency. In the spectral dimension of the excitation profile, harmonic excitation sidebands arise with a periodicity of 1/t. In this case, the radiofrequency pulse is applied only during the positive gradient lobes. If the radiofrequency pulse were applied during both the positive and negative gradient lobes, main harmonics would appear at period 2/t and secondary harmonics at 1/t [START_REF] Bernstein | CHAPTER 5 -SPATIAL RADIOFREQUENCY PULSES[END_REF], but this is difficult to achieve in a robust manner with high-resolution NMR hardware and is not considered here. Because of these excitation sidebands, care has to be taken when setting up the target profile to avoid the unintended excitation of harmonic resonance frequencies. As mentioned by Shrot and Frydman [START_REF] Shrot | Spatial/spectral encoding of the spin interactions in ultrafast multidimensional NMR[END_REF], the encoding gradient oscillation period can be slightly modified to make sidebands coincide with empty regions of the spectra. Nevertheless, the implementation of SPSP encoding can be difficult for crowded spectra with large spectral widths and limited empty regions.

Excitation sidebands are detrimental and the question is how to increase the "excitation bandwidth" accessible without the presence of harmonic excitation. Intuitively, a parallel can be done with the problem of the maximum detection spectral width accessible in the frequency dimension of the (k, W) domain that constitutes an UF2DNMR spectrum. It is conveniently solved by an interleaved acquisition where N successive scans are recorded to sample interleaved (k, t) trajectories. In practice, this interleaving is achieved by adding a preacquisition delay incremented by 2Ta/N. In the combined dataset, the effective dwell time is divided by a factor N, resulting in a multiplication of the spectral width by the same factor. By analogy, a multi-scan solution to suppress excitation sidebands would be to interleave N segments in excitation (k, t) domain. This is the solution explored in this paper, which is expected to work in the regime where the small-tip-angle approximation is valid, and relies on the linear relationship expressed in Eq. ( 8) [START_REF] Hardy | 31P Spectroscopic localization using pinwheel NMR excitation pulses[END_REF][START_REF] Panych | Selection of high-definition 2D virtual profiles with multiple RF pulse excitations along interleaved echo-planar k-space trajectories[END_REF].

Note that the properties of excitation sidebands can also be modified by adding a phase shift of p between subpulses during the design [START_REF] Schick | Simultaneous highly selective MR water and fat imaging using a simple new type of spectral-spatial excitation[END_REF][START_REF] Zur | Design of improved spectral-spatial pulses for routine clinical use[END_REF]. The resulting excitation profile appears shifted along the frequency axis and the central band is suppressed. This method is well adapted for selection/suppression in MRI but is not appropriate for UF2DNMR spectroscopy, where the central band is the targeted one.

Interleaved spatial/spectral encoding

Interleaved SPSP encoding may be described as the linear combination of data obtained from different segmented trajectories in excitation (k, t) domain. The combined grid has a finer time sampling and results in a larger excitation bandwidth. Interleaved SPSP encoding can also be described as a form of phase cycling, in which excitation sideband signals in the spectral dimension of the (k, W) domain are suppressed through the combination of N successive scans. When interleaved excitations are performed with a time offset Dt between shifted trajectories in excitation (k, t) domain, the phase difference between scans Df for a resonance frequency n in the rotating frame is:

Δ𝜙 = 2𝜋 𝜈 Δ𝑡. (9)
As excitation sidebands appear with a period 1/t when the radiofrequency pulse is applied only during the positive encoding gradient lobes (with t oscillation period), the phase difference for the jth frequency harmonic is:

Δ𝜙 = 2𝜋 𝑗 C r Δ𝑡. ( 10
)
When N interleaved excitations are used, trajectories in excitation (k, t) domain are shifted by a time offset Dt set to t/N to obtain:

Δ𝜙 = 2𝜋 𝑗/𝑁. ( 11
)
As a consequence, the principal band (j = 0) undergoes no phase difference between scans as well as all the jth band multiple of N. All other frequency harmonics undergo a phase cycling among Eq. [START_REF] Daniele | Single-Scan Multidimensional NMR Analysis of Mixtures at Sub-Millimolar Concentrations by using SABRE Hyperpolarization[END_REF]. At the end, the sum of all N scans leads to the cancellation of all but one out of N sidebands in the spectral dimension of the (k, W) domain. This multi-scan approach requires an adaptation of the pulse design to take into account time shifts used for interleaving. In practice, the N scans are recorded with N different SPSP pulses with interleaving delays before and after the pulse. The pulses are designed with the same target profile but with different interleaved (k, t) trajectories time shifted by multiple of t/N. For the nth scan, the delays di and di' respectively placed before and after the pulse are set according to: 𝛿 ) = (𝑛 -1) 𝜏/𝑁, [START_REF] Pelupessy | Improving resolution in single-scan 2D spectroscopy[END_REF] and

𝛿 ) l = (𝑁 -𝑛) 𝜏/𝑁. ( 13
)
This principle is presented in Fig. 3c with two interleaved SPSP excitations. In this particular case, the time offset of the (k, t) trajectories between the two scans is t/2. Two different SPSP pulses are designed to match the two interleaved (k, t) trajectories in Fig. 3e. Here, the phase difference of the jth sideband between the two scans is Δ𝜙 = 𝑗𝜋. In the simulated response profile of transverse spin magnetization |MXY(z, W)| obtained after the sum of the two scans (see Fig. 3g), the signals of the principal band (j = 0) and of all jth bands multiple of 2 interfere in a constructive manner. All other jth sideband signals undergo a phase cycling with p steps and interfere in a destructive manner. The sum of the two scans leads to the cancellation of one out of two harmonics in the spectral dimension of the (k, W) domain.

In applying this principle with the appropriate number of interleaved scans, a selective excitation of signals over a large spectral width can be performed with an efficient removal of periodic excitation sidebands. The pulse sequence in Fig. 4a shows the implementation of the interleaved SPSP encoding. The role of this basic pulse sequence is to perform a single excitation without mixing to yield a 2D spectrum with only diagonal peaks. Since the SPSP pulse is applied directly at the beginning of the pulse sequence, the interleaving delay di before the SPSP pulse is unnecessary. To put interleaved SPSP encoding into practice, one just has to replace the classical encoding scheme used in UF2DNMR experiments by a series of custom designed SPSP pulses, with associated interleaving delays. The use of a p/2 angle for SPSP pulses is preferred to a p angle. In the latter case, the small-tip-angle regime is no longer valid, leading to a degradation of the encoding efficiency. In theory, any UF2DNMR sequence can be modified with this logic. 13). The narrow and wide black rectangles stand for hard p/2 and p pulses. Black shaped pulses represent the RF part of p/2 SPSP pulses. The associated SPSP encoding gradient strength Ge, oscillation period t and number of cycles NCycle are provided in captions of Figs.

5-7 as well as EPSI acquisition parameters (acquisition gradient strength Ga and duration Ta, and number of repetitions NLoop).

Results and discussion

Interleaved SPSP excitation was first exemplified with cyclosporin A whose 1 H spectrum, shown in Fig. 5a, has a large frequency dispersion over 8 ppm. An SPSP pulse was designed to select and encode six resonances distributed along the spectrum at six different positions on the k-axis. The pulse sequence in Fig. 4a was used to perform a single-scan mixingless correlation experiment. Fig. 5c shows the resulting cross-peak-free spectrum where the targeted signals appear at targeted k-axis indexes in the spatially encoded dimension, together with signals whose frequency coincides with excitation sidebands of the SPSP pulse. Weak signals also arise from the intense methyl singlets that are partially excited even if they lie outside of the targeted spectral regions.

For this experiment, an encoding gradient with 1 ms oscillation period was used, yielding excitation sidebands every 1000 Hz. Gradient hardware instabilities did not permit to reduce the period under a value of 1 ms. In these conditions, to cover the whole 4000 Hz spectral width, a total of four interleaved scans were necessary to eliminate unwanted excitation sidebands. The sum of all scans yielded the spectrum in Fig. 5d, totally free of sidebands, with only signals targeted by the SPSP encoding with (k-space position/offset) pairs. In the four different interleaved scans, sidebands were phase cycled as shown in Fig. 5b for the signal at k-axis index 1 in the spatially encoded dimension. This signal close to 8 ppm corresponding to the principal band kept a constant phase in the four scans. The first frequency harmonic coincided with an empty region of the spectrum close to 6 ppm. The second one close to 4 ppm was phase cycled with p steps and the third one close to 2 ppm was phase cycled with 3p/2 steps as expected theoretically. All k-axis index traces extracted from the 2D spectra presented in Fig. 5 are provided as supplementary information (Fig. S1). Note that artefactual signals were observed experimentally at k = 0 in the spatially encoded dimension. For this reason, targeted peaks were all assigned to positive values of k and rephasing gradient following SPSP pulse and prephasing gradient before acquisition were omitted. After the interleaved excitation was validated, the principle was extended to a selective UF TOCSY sequence (Fig. 4b). The goal of this pulse sequence was to transfer, by isotropic mixing, the spatially encoded spin magnetizations of all the nuclei initially excited by the SPSP pulse within the entire spin systems to which they belong. The TOCSY experiment was a particularly relevant case study to test the performance of the interleaving excitation scheme because a large number of signals resulted from the use of isotropic mixing. The pulse sequence in Fig. 4b yielded the spectrum in Fig. 6a within a single scan using the previously defined SPSP encoding conditions. For each k-axis index in the spatially encoded dimension, a mix of the selective TOCSY spectrum from signal of interest and selective TOCSY spectra from all excitation sidebands was observed. The sum of the four interleaved scans yielded the spectrum in Fig. 6b, containing only the individual sub-spectra of six residues of cyclosporin A. The result was equivalent to a set of selective 1D-TOCSY experiments [START_REF] Kessler | Multi-dimensional NMR experiments using selective pulses[END_REF], but with the UF2DNMR SPSP encoded approach all experiments were recorded at the same time. All k-axis index traces extracted from Fig. 6b are provided as supplementary information (Fig. S2) with full assignment of 1 H resonances. The best alternative selective strategy in terms of signal-tonoise ratio is Hadamard spectroscopy [START_REF] Kupče | Hadamard NMR spectroscopy[END_REF] which is based on the multiplexed recording of selective experiments [START_REF] Bircher | Improved sensitivity for the 2D analogue of the 3D HOHAHA-COSY experiment[END_REF][START_REF] Blechta | Multi-site Hadamard NMR spectroscopy[END_REF][START_REF] Schraml | Hadamard 1D 1H TOCSY and its application to oligosaccharides[END_REF][START_REF] Plainchont | An alternative scheme for the multiplexed acquisition of 1D and 2D NMR spectra[END_REF]. After the demonstration of the use of interleaved SPSP encoding for 1 H-1 H correlation experiments, it was exemplified in a heteronuclear UF HSQC experiment with a mixture of metabolites. In the pulse sequence in Fig. 4c, the SPSP pulse was applied on the 13 C nuclei and substituted the p/2 pulse responsible for the coherence transfer in the INEPT block. The interleaving delay di before the SPSP pulse was not necessary in the absence of transversal magnetization before the pulse. An SPSP pulse was designed to select the eight signals of the 13 C spectrum at seven different positions in the k-axis. The two C6 13 C resonances of a and b forms of D-glucose were not distinguished due to poor separation (20 Hz). The pulse sequence in Fig. 4c yielded the spectrum in Fig. 7a within a single scan. Each 1 H resonance correlated with multiple carbon k-axis indexes. The reason is that the encoding gradient oscillation period was 1 ms, thus excitation sidebands occurred every 1000 Hz and coincided with almost every 13 C resonances. The interleaved SPSP encoding of 13 C resonances required a total of eight scans to cover the frequency dispersion of 8000 Hz with complete suppression of unwanted excitation harmonics. The sum of all scans yielded the spectrum in Fig. 7b, with only one 13 C resonance assigned to each k-axis index according to the (k-space position/offset) pairs used to design the SPSP pulse, resulting in one HSQC correlation per k-axis index. The SPSP encoding of carbon nuclei required more gradient power because of the lower g value. In this example, 78% of the maximum gradient strength was used to encode seven offsets. It follows that this method has a limitation in the number of resonances to encode while keeping a sufficient resolution of signals in the spatially encoded dimension.

Conclusions

Spatial/spectral encoding is a potential solution to address the spectral width and resolution limitations that are inherent to the linear encoding approach typically used in UF2DNMR experiments. Single-scan 2D spectra obtained using SPSP encoding may however display artefactual signals because of the excitation sidebands of SPSP pulses. We have carried out a theoretical and numerical analysis to develop a method to efficiently eliminate unwanted sidebands. We propose a multi-scan interleaving strategy in the excitation (k, t) domain, using a set of custom designed SPSP pulses. The proposed scheme is illustrated with a UF selective TOCSY strategy to assign the 1 H NMR resonances of some individual residues of the cyclosporin A peptide, and a UF HSQC spectrum of a metabolite mixture. The interleaved SPSP encoding procedure can benefit to the analysis of samples with large chemical shift dispersion, the use of UF2DNMR at high magnetic fields and the structural analysis of individual components of small molecule mixtures.

Experimental

All numerical simulations were performed with the Spinach package [START_REF] Hogben | Spinach -A software library for simulation of spin dynamics in large spin systems[END_REF][START_REF] Kuprov | Large-scale NMR simulations in liquid state: A tutorial[END_REF] for MATLAB. The magnetic field was set to 11.7 T. The sample dimension was fixed to 15 mm in the spatial dimension using a linear grid of 500 points. A spectral region of 15652 Hz was considered in the spectral dimension, using a linear grid of 500 points (in practice, this was implemented using a second spatial dimension, with a length of 4 mm and a constant gradient of 0.091 T.m - 1 ). The SPSP encoding pulse was designed with the following parameters: Ge = 0.195 T.m -1 , t = 600 µs and NCycle = 25. Each radiofrequency subpulse was defined with 100 points and applied during only the positive gradient lobes.

The cyclosporin A sample was prepared by dissolution in C6D6 at a 25 mM concentration. The metabolite mixture was prepared with [2- 13 All spectra were recorded at 303 K on a 500 MHz Avance III Bruker spectrometer equipped with a cryogenic 1 H/ 13 C probe with a z-axis gradient coil. UF pulse sequence details are given in captions of Figs. 4567. All gradient values are expressed as percents of the hardware maximum intensity (65 G.cm -1 ). SPSP pulses were designed with a MATLAB script performing a non-uniform Fourier transform from a target pattern with an algorithm based on conjugate gradients which converged within 3 iterations [START_REF] Shrot | Spatial/spectral encoding of the spin interactions in ultrafast multidimensional NMR[END_REF]. Radiofrequency and gradient pulses were exported in the Bruker file format.

The processing of UF spectra was performed with a MATLAB script. A spatial Gaussian apodization was applied in the spatially encoded dimension [START_REF] Giraudeau | Sensitivity and lineshape improvement in ultrafast 2D NMR by optimized apodization in the spatially encoded dimension[END_REF]. A sinusoidal apodization was applied in the frequency dimension. Zero-filling was applied in both dimensions.

Fig. 1 .

 1 Fig. 1. Spatial encoding schemes in ultrafast 2D NMR and spatial winding of the magnetization along the sample, with the corresponding echo position in the k-space dimension. After spatial encoding, each of the three chemical-shift offsets is associated to a specific helix pitch. The unwinding of each helix during the acquisition gradient Ga leads to an echo when the spatially dependent phase of the transverse spin magnetization is refocused along the sample. (a) Linear encoding with a pair of adiabatic frequency-swept pulses applied during a bipolar gradient produces a linear dependence of k-space position versus chemical shift value (b). (c) The signal consisting in a series of echoes is decoded in the k-space range during the acquisition gradient Ga. (d) Arbitrary encoding with spatial/spectral pulses allows arbitrary position in the k-space for each chemical shift (e), and an optimally reduced kspace acquisition range with equally spaced echoes (f).

Fig. 3 .

 3 Fig. 3. (a) Amplitude of the target (z, W) domain profile 𝑀 TU V-WXEY with only one on-resonance selection in the spectral dimension. (f, g) Amplitude of the simulated response of transverse spin magnetization to SPSP encoding with gradient oscillation period t = 600 µs and radiofrequency pulse applied during only the positive gradient lobes. (left) After the application of a single SPSP excitation (b) corresponding to a given k-space trajectory (d), excitation sidebands occur every 1667 Hz (1/t) on the simulated |MXY(z, W)| profile (f). (right) With two interleaved scans (c), each SPSP pulse is designed to match two interleaved k-space trajectories with a time offset equal to t/2 (e), and the sum of the two scans leads to the elimination of half of the excitation sidebands in the simulated |MXY(z, W)| profile (g). With two interleaved SPSP excitations, the apparent excitation sideband periodicity is doubled (every 3333 Hz in this example).

Fig. 4 .

 4 Fig. 4. Pulse sequences for ultrafast 2D with interleaved SPSP encoding. (a) Basic ultrafast mixingless sequence. (b) Ultrafast TOCSY sequence. (c) Ultrafast HSQC sequence. The variable interleaving delay di' is set according to Eq. (13). The narrow and wide black rectangles stand for hard p/2 and p

Fig. 5 .

 5 Fig. 5. (a) Structure of cyclosporin A and its 1 H NMR spectrum. The arrows mark the six 1 H resonances selected to design the SPSP pulse. (b) Traces extracted at k-axis index 1 (selection of NH signal of Abu-2 residue) from the individual scans and from the sum of four interleaved scans before magnitude treatment. (c) Single-scan SPSP encoded ultrafast selective spectrum of cyclosporin A with selection of six 1 H resonances, without sideband suppression. (d) SPSP encoded ultrafast selective spectrum of cyclosporin A obtained with the sum of four interleaved scans with selection of six 1 H resonances. The ultrafast spectra were obtained with sequence in Fig. 4a. Ge = 16, t = 1 ms, NCycle = 35. Ga = 40, Ta = 223.2 µs, NLoop = 128, with 2 interleaved acquisitions.

Fig. 6 .

 6 Fig. 6. (a) Single-scan SPSP encoded ultrafast selective TOCSY spectrum of cyclosporin A with selection of six 1 H resonances, without sideband suppression. (b) SPSP encoded ultrafast selective TOCSY spectrum of cyclosporin A obtained with the sum of four interleaved scans with selection of six 1 H resonances. The spectra were obtained with the pulse sequence in Fig. 4b and the same signals of interest as in Fig. 5 were selected by SPSP encoding. MLEV-16 duration is 70 ms with 6 kHz intensity. Ge = 16, t = 1 ms, NCycle = 35. Ga = 40, Ta = 223.2 µs, NLoop = 128, with 2 interleaved acquisitions. G1/G2 = 80/-80 (1 ms).

Fig. 7 .

 7 Fig. 7. (a) Single-scan SPSP encoded ultrafast HSQC spectrum of a metabolite mixture with selection of seven 13 C resonances, without sideband suppression. (b) SPSP encoded ultrafast HSQC spectrum of the same metabolite mixture obtained with the sum of eight interleaved scans with selection of seven 13 C resonances. The spectra were obtained with sequence in Fig. 4c. Ge = 78, t = 1 ms, NCycle = 25. Ga = 30, Ta = 389.6 µs, NLoop = 128. G1/G2/G3/G4 = 70/80/70/-20 (1 ms).

  C] Glycine, [2-13 C] L-alanine, [1-13 C] D-glucose, [2-13 C] D-glucose and [6-13 C] D-glucose dissolved in D2O at a 200 mM concentration.

Acknowledgments

This work was supported by the French National Research Agency (grant ANR-17-ERC2-0011), the French CNRS, the Université de Nantes, the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 801774), and the Région Pays de la Loire (Connect Talent).