Francis Comets
email: comets@math.univ-paris-diderot.fr

Michael Cranston

Overlaps and Pathwise Localization in the Anderson Polymer Model

Keywords: Brownian polymer, overlap, Malliavin calculus, parabolic Anderson model. AMS Subject classification numbers: Primary 60K35, 60K37; Secondary 60H05 60J65 82B44

We consider large time behavior of typical paths under the Anderson polymer measure. If P x κ is the measure induced by rate κ, simple, symmetric random walk on Z d started at x, this measure is defined as

where

) the normalizing constant. We establish that the polymer measure gives a macroscopic mass to a small neighborhood of a typical path as T → ∞, for parameter values outside the perturbative regime of the random walk, giving a pathwise approach to polymer localization, in contrast with existing results. The localization becomes complete as β 2 κ → ∞ in the sense that the mass grows to 1. The proof makes use of the overlap between two independent samples drawn under the Gibbs measure µ x κ,β,T , which can be estimated by the integration by parts formula for the Gaussian environment. Conditioning this measure on the number of jumps, we obtain a canonical measure which already shows scaling properties, thermodynamic limits, and decoupling of the parameters.

Introduction

In this paper we consider a polymer model related to the parabolic Anderson equation. In particular, we give quantitative bounds on the overlap of the polymer measure in terms of an inverse temperature parameter. This gives a quantitative expression for the extent to which the polymer measure concentrates its weight near a particular path at low temperature. The Anderson polymer model is a measure on simple, symmetric, continuous-time random walks influenced by a random field, W = {W x : x ∈ Z d }, of iid Brownian motions defined on some probability space (Ω, F t , Q) where the filtration is given by F t = σ({W x (s) : 0 ≤ s ≤ t, x ∈ Z d }). (Expectation with respect to Q will be denoted by E.) This measure and related quantities have motivated a huge number of research papers from many different perspectives. In the seminal reference [START_REF] Carmona | Parabolic Anderson problem and intermittency[END_REF] the model is used to give a mathematical account to intermittency, i.e. the existence of spots where most of the mass is concentrated. Those spots correspond to favorable configurations of the field. The model is also very appealing, with time-space iid environment (the Brownian increments) replaced by the configuration of an interacting particle system [START_REF] Gärtner | Intermittency in a catalytic random medium[END_REF], modeling a chemical reaction with moving catalysers. It has non trivial large deviation properties [START_REF] Cranston | On large deviations for the parabolic Anderson model[END_REF], as a particular random growth model. The one-dimensional totally asymmetric case, where the walker only jumps to the right, has a lower complexity than the symmetric case, as shown by the computations of annealed Lyapunov exponents [START_REF] Borodin | On moments of the parabolic Anderson model[END_REF]; In this case an explicit solution was given in [START_REF] Moriarty | On the free energy of a directed polymer in a Brownian environment[END_REF], with the strongly asymmetric case as a small perturbation [START_REF] Flores | Asymmetric directed polymers in random environments[END_REF]. We consider the symmetric case and, like [START_REF] Cadel | Sharp asymptotics for the partition function of some continuous-time directed polymers[END_REF], [START_REF] Carmona | Strong disorder implies strong localization for directed polymers in a random environment ALEA Lat[END_REF], [START_REF] Rovira | On the Brownian-directed polymer in a Gaussian random environment[END_REF] and [START_REF] Seppäläinen | Bounds for scaling exponents for a 1 + 1 dimensional directed polymer in a Brownian environment[END_REF], our focus is to view the parabolic Anderson model as a particular model of directed polymers in random medium, and a host of other references which may be found in [START_REF] Comets | Probabilistic analysis of directed polymers in a random environment: a review[END_REF] for a general picture.

In order to describe this model, start with the measure P x κ to be the measure on the canonical probability space, D([0, ∞), Z d ) of right continuous paths which possess left limits everywhere, with a finite number of jumps of size one only on any finite interval. These are the typical sample paths of the simple symmetric rate κ random walk. Here, as is usual, P x κ (X 0 = x) = 1 and with respect to P x κ , the canonical process X(t, ω) = ω(t), ω ∈ D ∞ = D([0, ∞), Z d ) is a Markov process with infinitesimal generator κ∆ where ∆ is the discrete Laplacian defined by ∆f (x) = 1 2d ||y-x||=1 (f (y)f (x)). The Anderson polymer model is the Gibbs measure on D T = D([0, T ], Z d ) defined by

µ x κ,β,T (f ) = Z κ,β,T (x) -1 E x κ f exp β T 0 dW X(s) (s) .
for bounded measurable f : D T → R. The model has three parameters, the inverse temperature β ∈ R measuring the fluctuations of the environment, the polymer length T and the diffusivity κ ∈ (0, ∞) of the path under the a priori measure P x κ . The partition function Z κ,β,T (x) is given by

Z κ,β,T (x) = E x κ exp β T 0 dW X(s) (s) .
By the Feynman-Kac formula,

u(t, x) = E x κ exp β t 0 dW X(t-s) (s)
is the solution of the time-dependent parabolic Anderson equation (or stochastic heat equation)

u(t, x) = 1 + κ t 0 ∆u(s, x)ds + β t 0 u(s, x) • dW x (s), (1)
where •dW denotes the Stratonovich differential of W . In the space-continuous model, the logarithm of u formally solves the Kardar-Parisi-Zhang equation (KPZ), which is expected to be the scaling limit of our model in a weak noise limit. In dimension d = 1, its distribution has been recently computed in [START_REF] Amir | Quastel: Probability Distribution of the Free Energy of the Continuum Directed Random Polymer in 1+1 dimensions[END_REF], and a better understanding of the KPZ universality class is being achieved; see [START_REF] Corwin | The Kardar-Parisi-Zhang equation and universality class[END_REF] for a review. The functions u(t, x) and Z κ,β,t (x) thus have the same distribution and we will make use of the properties of u(t, x) derived in [START_REF] Cranston | Lyapunov exponents for the parabolic Anderson model[END_REF] and apply them to Z κ,β,t (x). We shall use the notation E κ , µ κ,β,T and Z κ,β,T when x = 0. By spatial homogeneity of the field W we may confine ourselves to the study of the case x = 0. Our results concern the behavior of µ κ,β,T for β 2 /κ and T large in arbitrary dimension. It will be shown that this measure concentrates its mass near a particular favorite path as β 2 /κ → ∞. In a sense this means that there is a channel in the media in which most of the polymer paths reside. The establishment of the concentration is done by examining the overlap defined as

J κ,β,T ≡ 1 T T 0 µ ⊗2 κ,β,T (X(t) = X(t))dt (2)
where µ ⊗2 κ,β,T is the product measure of µ κ,β,T with itself. The paths X and X appearing in (2) are thus independent samples of paths drawn according to the measure µ κ,β,T . The quantity J κ,β,T is, for two independent samples sharing the same environment, the proportion of time spent together. We also study another version of the overlap. Define

I κ,β,T = 1 T T 0 µ ⊗2 κ,β,t (X(t) = X(t))dt . (3)
This version of the overlap measures the amount of time up to T that the endpoint of independent samples drawn with respect to the measure µ κ,β,t agree.

In statistical mechanics, counterparts of these overlaps can be found, for the Sherrington-Kirkpatrick model and other ones for disordered systems. Coming via integration by parts, the first overlap revealed most successful in the last decade ( [START_REF] Guerra | Sum rules for the free energy in the mean field spin glass model[END_REF], [START_REF] Guerra | The thermodynamic limit in mean field spin glass models[END_REF], [START_REF] Talagrand | Spin Glasses: A Challenge for Mathematicians. Cavity and Mean Field Models[END_REF], [START_REF] Talagrand | The Parisi formula[END_REF]) to study the low temperature regime.

Below, we shall see that the overlap in our polymer model arises in a similar fashion. By taking the logarithmic Malliavin derivative of the partition function with respect to W x (t) one arrives at J κ,β,T . By taking the logarithmic Itô derivative of the partition function one arrives at I κ,β,T . Both these overlaps take values in [0,1], they vanish in the limit of large T 's under the free measure, i.e. for β = 0, and positivity of each one implies a specific form of localization. It is known that positivity of I κ,β,T is equivalent to positivity of the difference between the annealed and the quenched free energies. Here, we show that positivity of J κ,β,T essentially amounts to this difference being strictly increasing as a function of |β|. Then, using the logarithmic asymptotics of the partition function from [START_REF] Cranston | Lyapunov exponents for the parabolic Anderson model[END_REF], we find that both versions of the overlap converge to 1 as β 2 /κ → ∞, hence achieving its upper bound. This allows us to give a precise account of path localization by identifying the favourite end-point, and the favourite path for the polymer. We introduce here a sequence of measurable functions y * T : [0, T ] → Z d such that the proportion of time when the polymer is equal to the favourite path is positive in the localized phase, and converges to 1 as

β 2 /κ → ∞: lim inf T →∞ Eµ κ,β,T 1 T T 0 δ 0 X(t) -y * T (t) dt -→ 1 as β 2 /κ → ∞. (4)
The function y * T , defined in [START_REF] Marquez-Carreras | A model of continuous time polymer on the lattice[END_REF], depends on κ, β and on the environment, it is called the favourite path. Our statement improves on the literature on polymers by concerning the path itself, not only the terminal location at time T as in [START_REF] Carmona | On the partition function of a directed polymer in a random environment[END_REF] and [START_REF] Comets | Directed polymers in a random environment: path localization and strong disorder[END_REF]. Also, we can obtain, in this model, complete localization in the limit, in the sense of (4) where the right hand side is equal to its maximum value 1.

We now mention related results on complete localization. Strong concentration for the directed polymer in a random environment for parabolic Anderson model (space dependent only) with a Pareto potential was established in [START_REF] König | A two cities theorem for the parabolic Anderson model[END_REF]. The main difference is that there, the favourable sites in the environment have a simple characterization in terms of the potential. In the time-discrete case with heavy tailed potentials (timespace dependent), see [START_REF] Auffinger | Directed polymers in random environment with heavy tails[END_REF] for similar conclusions. When the tails are less heavy, the favourite corridors can no longer be characterized by maxima of the potential, they depend in a much more subtle manner on the environment; Though they are not anymore explicit, site localization can still be proved [START_REF] Vargas | Strong localization and macroscopic atoms for directed polymers[END_REF]. Note that in the discrete case, only little is known on the random geodesics [START_REF] Newman | Topics in disordered systems[END_REF] in first passage percolation, which are the zero-temperature favourite paths. In our parabolic Anderson model, the potential has strong decay, but we can prove strong localization. For the solution of the KPZ equation in one dimension, the distribution of the favourite end-point has been recently computed in [START_REF] Flores | Remenik Endpoint distribution of directed polymers in 1+1 dimensions[END_REF], it is the arg max of an Airy 2 process minus a parabola.

The parabolic Anderson model, compared to other directed polymers, has some nice scaling properties, it also decouples some parameters due to its Poissonian structure. With a Gaussian potential, it allows a simple integration by parts formula, which comes from the Malliavin calculus.

The organization of this paper is as follows, in Section 2 we present our main results. Section 3 contains preliminaries from the Malliavin calculus giving a first version of the overlap. Section 4 gives an Ito calculus derivation of the second version of the overlap, and contains the proofs for localization. In Section 5 we prove the results for the distribution of the number of jumps under the polymer measure, and regularity properties of the favourite end-point and path in Section 6.

Main Results

From now on, we adopt the notation

H T (γ) = T 0 dW γ(s) (s), γ ∈ D T .
(5)

Thermodynamic limits

Our first goal is to determine the large deviation rate for the number of jumps of paths with respect to the measure µ κ,β,T . We begin with the exponential growth rate of Z κ,β,T conditioned on the number of jumps of the process X up to time T which we denote by N(T, X).

Proposition 2.1. For r ≥ 0, the limit

Γ(β, r) = lim T →∞ T -1 ln E κ [exp{βH T (X)}|N(T, X) = [rT ]] (6)
exists a.s. and in L p , p ∈ [1, ∞). The limit is deterministic, symmetric and convex in β, continuous in (β, r), independent of κ, and satisfies the scaling relation

Γ(β, r) = aΓ(a -1/2 β, a -1 r), a > 0. (7)
Let I κ be the Cramér transform (i.e. the large deviation rate function) of the Poisson distribution with parameter κ,

I κ (r) = r ln(r/κ) -r + κ, r ≥ 0. ( 8
)
We will see that the function r → Γ(β, r) -I κ (r) is concave on R + , and tends to -∞ as r → +∞.

Proposition 2.2 (Free energy). The limit

Ψ(κ, β) = lim T →∞
T -1 ln Z κ,β,T exists a.s. and in L p , p ∈ [1, ∞), and is equal to

Ψ(κ, β) = sup{Γ(β, r) -I κ (r); r ≥ 0}. (9)
In particular, for all positive a,

Ψ(κ, β) = aΨ(a -1 κ, a -1/2 β) = β 2 Ψ(β -2 κ, 1). = κΨ(1, κ -1/2 β). (10)
From this we obtain the quenched large deviation rate function for the distribution of the number of jumps of the polymer.

Theorem 1 (Large deviations). Define I κ,β to be the convex function

I κ,β (r) = -Γ(β, r) + I κ (r) + Ψ(κ, β). Then lim T →∞,n/T →r T -1 ln µ κ,β,T N(T, X) = n = -I κ,β (r), a.s.. (11)
Moreover, for a.e. realization of the environment W = {W x (•) : x ∈ Z d }, and all subsets B ⊂ R + ,

-inf B o I κ,β (r) ≤ lim inf T →∞ T -1 ln µ κ,β,T N(T, X)/T ∈ B ≤ lim sup T →∞ T -1 ln µ κ,β,T N(T, X)/T ∈ B ≤ -inf B I κ,β (r). Remark 2.1. (i) Annealed bound. By Jensen's inequality, it is readily checked that both Ψ(κ, β) ≤ β 2 /2 and Γ(β, r) ≤ β 2 /2
hold. These, in addition to [START_REF] Carmona | Parabolic Anderson problem and intermittency[END_REF], imply that

Ψ(κ, β) = β 2 /2 ⇐⇒ Γ(β, κ) = β 2 /2, (12)
and, in such a case, I κ,β (r) has a unique minimum at r = κ.

(ii) Weak versus strong disorder. From ( 31), [START_REF] Vargas | Strong localization and macroscopic atoms for directed polymers[END_REF] it follows that

Ψ(κ, β) = β 2 /2 ⇐⇒ β 2 /κ ≤ Υ c , (13)
for some critical value Υ c ∈ [0, ∞) depending only on the dimension; finiteness of Υ c can be seen from ( 25) below.

In dimension d ≥ 3, it is known by second moment method [START_REF] Carmona | Strong disorder implies strong localization for directed polymers in a random environment ALEA Lat[END_REF] that Z κ,β,T exp{-T β 2 /2} converges to a positive limit, so that the equality holds in the left member of ( 12) for β 2 /κ small. Hence, Υ c > 0 in that case.

In dimension d = 1, 2, it has recently been proved [START_REF] Bertin | Very strong disorder for the parabolic Anderson model in low dimensions[END_REF] that Υ c = 0, by extending the techniques introduced in [START_REF] Lacoin | New bounds for the free energy of directed polymers in dimension 1+1 and 1+2[END_REF] for discrete models.

Overlaps and phase transition

Recall the definitions (3) of I κ,β,T and (2) of J κ,β,T . Even though the quantity J κ,β,T places more restriction on the paths, measuring the fraction of time together from 0 to T with respect to µ ⊗2 κ,β,T , for large β 2 /κ it is essentially the same size as I κ,β,T . The advantage of J κ,β,T is that it involves a single Gibbs measure and therefore contains pathwise information. We will prove the following result in Sections 3 and 4. We remark that (15) was also proven in [START_REF] Carmona | Strong disorder implies strong localization for directed polymers in a random environment ALEA Lat[END_REF]. Also, discrete time versions of (17) were established in [START_REF] Carmona | On the partition function of a directed polymer in a random environment[END_REF] and in [START_REF] Lacoin | New bounds for the free energy of directed polymers in dimension 1+1 and 1+2[END_REF].

Proposition 2.3. (Overlaps) (i) For all β and κ the limit

I κ,β,∞ = lim T →∞ I κ,β,T (14)
exists almost surely and is nonrandom, and is equal to

I κ,β,∞ = 1 - 2 β 2 Ψ(κ, β). (15)
(ii) The limit

J κ,β,∞ = lim T →∞ E [J κ,β,T ] (16)
exists for all κ and all β except for an at most countable set of values of β 2 /κ, and

J κ,β,∞ = 1 -β -1 ∂ ∂β Ψ(κ, β). (17)
The first step for (i) is an Itô calculation, and for (ii) it's an integration by parts formula, which relies on the Malliavin calculus.

We say that a function

f : R → R is non-constant around b ∈ R if, in all neigh- borhoods of b there is some b ′ with f (b ′ ) = f (b). Define the subsets of R, NC 1 = b : β 2 2 -Ψ(1, β) is non-constant around b , NC κ = κ 1/2 NC 1 .
By [START_REF] Comets | Directed polymers in a random environment: path localization and strong disorder[END_REF], β → β 2 2 -Ψ(κ, β) is non-constant around every point in NC κ , and only there. Note that, by [START_REF] Lacoin | New bounds for the free energy of directed polymers in dimension 1+1 and 1+2[END_REF], this function is non-decreasing as |β| is increased. Also, if β is not in the subdifferential of Ψ(κ, •) at β, then β ∈ NC κ .

Corollary 2.1. (i) We have, with Υ c from (13),

I κ,β,∞ > 0 ⇐⇒ β 2 /κ > Υ c .
(ii) Similarly, for all κ > 0,

Υ c = inf{β 2 /κ : J κ,β,∞ > 0}.
Moreover, if β ∈ NC κ , there exists a sequence β n converging to β such that J κ,βn,∞ > 0. Finally, if the derivative of Ψ(κ, •) at β exists and is different from β, then J κ,β,∞ > 0.

We also have asymptotic estimates on the overlaps.

Corollary 2.2. As β 2 /κ → ∞, we have (i) I κ,β,∞ = 1 - α 2 2 ln(β 2 /κ) (1 + o(1)),
and (ii)

J κ,β,∞ = 1 -O 1 ln(β 2 /κ) . (18)

End-point and path localization

For fixed κ, β we define the favourite end-point x * (t) for the polymer at time t by

x * (t) = arg max E κ [exp{βH t (X)}δ x (X(t))] : x ∈ Z d , (19)
taking the first argument x in the lexicographic order in case of multiple maxima. By definition, x * (t) maximizes the distribution µ κ,β,t (X(t) = x) of the end-point location.

Noting that for t ≤ T,

E κ [e βH T (X) ] = x∈Z d E κ [e βH T (X) δ x (X(t))]
is almost surely finite, we see that there is at least one maximizer. We will see in Proposition 2.4 (i) that, in fact, the spatial rate of decay is superexponential. We introduce for each T ≥ 0 a new object, the favourite path with time horizon T , defined as

y * T (t) = arg max E κ (exp{βH T (X)}δ x (X(t))) ; x ∈ Z d . (20)
At each time t ∈ [0, T ], it maximizes µ κ,β,T (X(t) = x); it does exist by a similar argument. The path t → y * T (t) is a.s. piecewise constant on [0, T ]. Indeed, if at time t there is only one maximizer x of f (x, t, T ) = E κ [exp{βH t (X)}δ x (X(t))] , then by continuity x will the be the unique maximizer for times close to t. In Proposition (2.4), (iv) below, we show that t → f (x, t, T ) is C 1 a.s. (denote by f' the derivative in t.) Now, assume that for x = y, f ′ (x, t, T )f ′ (y, t, T ) has a density. Then it is different from 0 a.s., and so one maximizer will win over the other ones for times close to t. We won't go into the details of the existence of the density, since we won't explicitly use the just stated claims. Otherwise, we know very little about the path t → y * T (t). We start with technical results. Denote

f (x, t) = E κ [exp{βH t (X)}δ x (X(t))] and f (x, t, T ) = E κ [exp{βH T (X)}δ x (X(t))] for short notations. Proposition 2.4. (i) For a < ln 2 there exists C 0 = C 0 (t, W) such that f (x, t) ≤ C 0 exp -{a|x| ln |x|}, x ∈ Z d .
(ii) For a < ln 2 there exists

C 1 = C 1 (T, W) such that f (x, t, T ) ≤ C 1 exp -{a|x| ln |x|}, x ∈ Z d , t ≤ T.
(iii) The function t → f (x, t) is almost surely Hölder continuous of every order less than 1/2.

(iv) The function t → f (x, t, T ) is almost surely of C 1 class on (0, T ).
We comment on some observations on the favourite attributes x * and y * T . Both depend on κ, β (also T for the second one) and on the environment W. Both have long jumps: x * / ∈ D ∞ a.s., and y * T / ∈ D T with positive probability for all T . By time continuity in (iii-iv) and since for a measurable F (ω, x), the (smallest) maximizer is a measurable function of ω, both (t, ω) → y * T (t) and (t, ω) → x * (t) are measurable functions from [0, T ]×Ω to Z d . Observe that the functions y * T and x * are equal at time t = T , but they are not related otherwise. In general,

x * (T /2) is different from y * T (T /2). The end- point process is F -adapted, with F = (F t ) t≥0 , F t = σ {W x (s) : 0 ≤ s ≤ t, x ∈ Z d } .
The other one , being only F T -measurable, is anticipating. Now, here is a fundamental difference between the two. The mapping t → x * (t) has oscillations at those times t when there are many maximizers: in view of Proposition 2.4 (iii), the set of jump times then looks locally like the set of zeros of Brownian motion. In contrast, from differentiability in Proposition 2.4 (iv), we see that t → y * T (t) has no oscillations. The favourite path is much smoother than the favourite end-point viewed as a process.

Coming to our main results, we prove that the polymer concentrates on the favorite end-point and path in the strong disorder region (loosely speaking), and overwhelmingly as β 2 /κ → ∞.

Theorem 2. (favourite site and path)

(i) If β 2 /κ > Υ c , there exists a constant C = C(β 2 /κ) > 0 such that lim inf T →∞ 1 T T 0 µ κ,β,t (X(t) = x * (t))dt ≥ C a.s. ( 21
)
On the contrary, for β 2 /κ ≤ Υ c , the left-hand side of ( 21) converges to 0. Finally, for all ǫ > 0, for all β 2 /κ large enough, we have a.s.,

(1 -ǫ) α 2 4 ln(β 2 /κ) ≤ lim inf T →∞ 1 T T 0 µ κ,β,t (X(t) = x * (t))dt ≤ lim sup T →∞ 1 T T 0 µ κ,β,t (X(t) = x * (t))dt ≤ (2 + ǫ) α 2 4 ln(β 2 /κ) . ( 22
)
(ii) For all κ > 0 and β ∈ NC κ , there exists a sequence

β n → β such that lim inf T →∞ Eµ κ,βn,T 1 T T 0 δ 0 (X(t) -y * T (t))dt > 0. ( 23
)
On the contrary, for

β 0 / ∈ NC κ , then lim T →∞ Eµ κ,β,T 1 T T 0 δ 0 (X(t) -y * T (t)
)dt = 0 for all β in a neighborhood of β 0 . Finally, for all ǫ > 0, for all β 2 /κ large enough, we have a.s.,

(1 -ǫ) α 2 4 ln(β 2 /κ) ≤ 1 -lim sup T →∞ Eµ κ,β,T 1 T T 0 δ 0 (X(t) -y * T (t))dt ≤ 1 -lim inf T →∞ Eµ κ,β,T 1 T T 0 δ 0 (X(t) -y * T (t))dt ≤ (2 + ǫ) α 2 4 ln(β 2 /κ) (24)
In words, the proportion of time the polymer sticks to the favourite path is positive in the long run if and only if the difference between annealed and quenched free energies is strictly increasing with |β|. Moreover, this proportion tends to 1 as the ratio β 2 /κ becomes large. In the case [START_REF] Moriarty | On the free energy of a directed polymer in a Brownian environment[END_REF] we say that path localization holds at (κ, β). In the case of [START_REF] Newman | Topics in disordered systems[END_REF], precisely when the overlap of the polymer with the favourite path tends to 1, we say that complete path localization takes place as

β 2 /κ → ∞.
Very little can be said so far on the favourite path y * T , which determines the corridor where most of the mass is concentrated. Both x * and y * T are complicated functions of the environment. However, from the above theorem, it is approximated, in the distance dist T (x, y) = T -1 T 0 δ 0 (x(t)-y(t))dt, within order of 1/ ln 2 (β 2 /κ)) accuracy, by a path γ from D T , with jump density T -1 N(T, γ) ∈ (r max (κ, β)ǫ, r max (κ, β) + ǫ) by Theorem 1. Since

κ -1 r max (κ, β) ∼ α 2 (β 2 /κ) 4 ln 2 (β 2 /κ) ,
the path y * T becomes wilder as β 2 /κ increases, but within certain limits. We end with a conjecture, which holds for polymer models on trees.

Conjecture 2.1.

NC κ = β : |β| ≥ (κΥ c ) 1/2 .

Asymptotics of Lyapunov exponents for parabolic Anderson model

The existence of the quenched Lyapunov exponent Ψ(κ, 1) is well-known in parabolic Anderson model literature, see for example [START_REF] Cadel | Sharp asymptotics for the partition function of some continuous-time directed polymers[END_REF], [START_REF] Carmona | Asymptotics for the almost sure Lyapunov exponent for the solution of the parabolic Anderson problem[END_REF], [START_REF] Carmona | Parabolic Anderson problem and intermittency[END_REF] and [START_REF] Cranston | Lyapunov exponents for the parabolic Anderson model[END_REF]. In [START_REF] Carmona | Asymptotics for the almost sure Lyapunov exponent for the solution of the parabolic Anderson problem[END_REF] and [START_REF] Cranston | Lyapunov exponents for the parabolic Anderson model[END_REF] it was shown that

Ψ(κ, 1) ∼ α 2 4 ln 1 κ , κ ց 0, ( 25
)
where ∼ means that the ratio of the two sides tends to one. A few words about the constant α are in order. Define the space of paths

Γ [0,n],n = γ ∈ D n : γ : [0, n] → Z d , N(γ, n) = n .
The superadditive functional

A n = sup γ∈Γ [0,n],n H n (γ)
is the supremum of a Gaussian field

H n (γ) : γ ∈ Γ [0,n],n
indexed by the set Γ [0,n],n . This set has a suitably bounded entropy, which allows the conclusion, see [START_REF] Carmona | Asymptotics for the almost sure Lyapunov exponent for the solution of the parabolic Anderson problem[END_REF] and [START_REF] Cranston | Lyapunov exponents for the parabolic Anderson model[END_REF], that there is a finite, positive constant α such that lim n→∞

1 n A n = α, a.s.
This is the constant α appearing in [START_REF] Nualart | The Malliavin calculus and related topics[END_REF]. Thus, by the scaling relation [START_REF] Comets | Directed polymers in a random environment: path localization and strong disorder[END_REF], it follows that

Ψ(κ, β) ∼ α 2 β 2 4 ln(β 2 /κ) , β 2 /κ → ∞ (26)
In particular, for k > 0,

lim β→∞ Ψ(β 2 e -kβ 2 , β) = α 2 4k .
We give a streamlined approach to [START_REF] Nualart | The Malliavin calculus and related topics[END_REF] and [START_REF] Rovira | On the Brownian-directed polymer in a Gaussian random environment[END_REF]. Obviously, since L p norms approach the L ∞ norm as p → ∞,

lim β→∞ 1 βT ln E κ exp{βH T (X)}|N(T, X) = [rT ] = sup γ∈Γ [0,T ],[rT ]
H T (γ)

≡ A T,r
By taking the limit as T → ∞ and interchanging the limits T → ∞ and β → ∞, a step to be justified later, see (50), we get

lim β→∞ β -1 Γ(β, r) = lim T →∞ T -1 A T,r . (27)
Now, by the Brownian scaling, A T,r L = √ rA T , making the previous limit equal to α √ r. Let r max (κ, β) be the set of maximizers of the right-hand side of ( 9). This set is a non-empty, compact interval included in R * + , and we conjecture it reduces to a single point. Observe by scaling that

r max (κ, β) = κr max (1, κ -1/2 β).
Parts (i) and (ii) of the following proposition were established in [START_REF] Carmona | Asymptotics for the almost sure Lyapunov exponent for the solution of the parabolic Anderson problem[END_REF], [START_REF] Cranston | Lyapunov exponents for the parabolic Anderson model[END_REF] and [START_REF] Marquez-Carreras | A model of continuous time polymer on the lattice[END_REF].

Ψ(κ, β) ∼ κ sup αβ r κ -I 1 (r) : r ≥ 0 ∼ α 2 β 2 4 ln(β 2 /κ) . (iii) As β 2 /κ → ∞, r max (κ, β) ∼ κ × α 2 (β 2 /κ) 4 ln 2 (β 2 /κ) .
Recall r max is defined as an interval, the last statement means that both endpoints are equivalent to the right-hand side. The behavior in (iii) for the typical number of jumps (r max (κ, β) >> κ) is drastically different from that in the weak disorder regime in Remark 2.1, (i).

Preliminaries from Malliavin calculus.

Express minus the Hamiltonian for a fixed path X as

H T (X) = x∈Z d T 0 dW x (t)δ x (X(t)).
This has the form

W (h) = x∈Z d T 0 h(t, x)dW x (t),
where h(t, x) depends on X by the relation h

(t, x) = δ x (X(t)). Obviously, h ∈ L 2 ([0, T ]×Z d ). The family {W (h) : h ∈ L 2 ([0, T ]×Z d )} is called a centered, isonormal
Gaussian family, and defines an abstract Wiener space as in [START_REF] Nualart | The Malliavin calculus and related topics[END_REF] or [START_REF] Üstünel | Transformation of measure on Wiener space[END_REF]. The Malliavin derivative DF of a square integrable random variable F defined on this space is, when it exists, a random element of L 2 ([0, T ] × Z d ), that we will view as a stochastic process DF = (D t,x F ) t,x indexed by time and space. The Malliavin derivative D t,x is heuristically equal to ∂ ∂(dWx(t)) and can be formally computed as such. The Malliavin derivative of H T (X) is thus the element of L 2 ([0, T ] × Z d ) defined by

D t,x H T (X) = δ x (X(t)).
Then taking f (y) = e y and applying the chain rule, we find the Malliavin derivative of f (βH T (X)) is given by

D t,x f (βH T (X)) = βf (βH T (X))δ x (X(t)).
Taking the average over paths X and then differentiating yields

D t,x Z κ,β,T = βE κ [δ x (X(t)) exp{βH T (X)}].
Note that we need to invoke not only linearity but also continuity to get this identity. Using again the chain rule, we obtain

D t,x ln Z κ,β,T = βµ κ,β,T (δ x (X(t))). ( 28
)
The crucial point is that

D ln Z κ,β,T 2
L 2 ([0,T ]×Z d ) = β 2 T 0 µ ⊗2 κ,β,T (X(t) = X(t))dt = β 2 T J κ,β,T ≤ β 2 T. (29)
The integration by parts formula

E[W (h)F ] = E x T 0 h D t,x F dt , (30)
is the central tool in Malliavin calculus, see [START_REF] Nualart | The Malliavin calculus and related topics[END_REF].

Lemma 3.1.

∂ ∂β E[ln Z κ,β,T ] = βT 1 -E[J κ,β,T ] . (31)
Proof. Differentiating inside the integral, we obtain

∂ ∂β E[ln Z κ,β,T ] = E[µ κ,β,T (H T (X))].
Then, we write

E [µ κ,β,T (H T (X))] =E µ κ,β,T x∈Z d T 0 dW x (t)δ x (X(t)) = x∈Z d EE κ e βH T (X) Z κ,β,T T 0 dW x (t)δ x (X(t)) = x∈Z d T 0 E κ E e βH T (X) Z κ,β,T dW x (t) δ x (X(t)) = x∈Z d T 0 E κ E D t,x e βH T (X) Z κ,β,T δ x (X(t)) dt
from Gaussian integration by parts, see [START_REF] Nualart | The Malliavin calculus and related topics[END_REF]. (A less pedestrian -though equivalent -computation is to apply directly the formula ( 30)). Then,

E [µ κ,β,T (H T (X))] = x∈Z d T 0 E κ E βD t,x H T (X) e βH T (X) Z κ,β,T - e βH T (X) Z κ,β,T E κ [D t,x H T ( X) e βH T ( X ] Z κ,β,T δ x (X(t)) dt = β x∈Z d T 0 E κ E δ x (X(t)) e βH T (X) Z κ,β,T - e βH T (X) Z κ,β,T E κ [δ x ( X(t))e βH T ( X) ] Z κ,β,T δ x (X(t)) dt = β x∈Z d T 0 E µ κ,β,T (δ x (X(t))) -µ κ,β,T (δ x (X(t))) 2 dt = βE T 0 1 -µ ⊗2 κ,β,T (X(t) = X(t)) dt = βT [1 -E[J κ,β,T ]] .
We will use the concentration of measure phenomenon in our analysis. The use of Malliavin calculus for concentration appeared in [START_REF] Rovira | On the Brownian-directed polymer in a Gaussian random environment[END_REF] in the study of polymers, and earlier in [START_REF] Toubol | High temperature regime for a multidimensional Sherrington-Kirkpatrick model of spin glass[END_REF] in the study of mean-field disordered systems.

Lemma 3.2 (Concentration).

Let A be a Borel subset of the path space with P κ x (A) > 0, and let

Z κ,β,T (A) = E κ [exp{βH T (X)}; A] .
Then, for all u > 0,

Q (| ln Z κ,β,T (A) -E ln Z κ,β,T (A)| ≥ u) ≤ 2 exp - u 2 2β 2 T .
Proof. Of course, Z κ,β,T (A) = Z κ,β,T when A is the full space. Following the above computations for Z κ,β,T , we see that the derivative is equal to

D t,x ln Z κ,β,T (A) = βµ κ,β,T (δ x (X(t))|A),
and, as in [START_REF] Talagrand | The Parisi formula[END_REF], its norm L 2 ([0, T ] × Z d ) is a.s. bounded by β 2 T . The lemma follows from Theorem B.8.1 in [START_REF] Üstünel | Transformation of measure on Wiener space[END_REF].

Next observe from (31) that

Ψ(κ, β) = lim T →∞ T -1 E [ln Z κ,β,T ] = lim T →∞ β 0 r [1 -E [J κ,r,T ]] dr. ( 32
) Define D ′ 1 = {β ∈ R : Ψ(1, β) is differentiable at β}, D 1 = D ′ 1 \ {0}, and
D κ = κ 1/2 D 1 .
By [START_REF] Comets | Directed polymers in a random environment: path localization and strong disorder[END_REF], β → Ψ(κ, β) is differentiable on D κ for all κ, and by convexity, and the complement of D κ is at most countable. By standard convexity arguments, the derivative

β(1 -E [J κ,β,T ]) of T -1 E [ln Z κ,β,T ] converges on D κ to (∂/∂β)Ψ(κ, β). Hence, the limit J κ,β,∞ := lim T →∞ E [J κ,β,T ] = 1 -β -1 ∂ ∂β Ψ(κ, β) (33)
exists for all κ and β ∈ D κ , with

Ψ(κ, β) = β 0 r 1 -J κ,r,∞ dr (34)
and moreover

β → β 1 -J κ,β,∞ is non decreasing. ( 35
) Proposition 3.1. As β 2 /κ → ∞, α 2 4 ln(β 2 /κ) (1 + o(1)) = 2 β 2 β 0 r 1 -J κ,r,∞ dr,
and

J κ,β,∞ = 1 -O 1 ln(β 2 /κ) .
Proof. The first statement is straightforward from ( 26) and (34). By (35), we have for β ≥ 0,

β 0 r 1 -J κ,r,∞ dr ≥ β β 2 r 1 -J κ,r,∞ dr ≥ β 2 4 1 -J κ,β/2,∞ .
Then, combining this with the first statement, we obtain

J κ,β/2,∞ ≥ 1 - α 2 2 ln β 2 κ (1 + o(1)),
yielding the same bound for J κ,β,∞ . This completes the proof.

An Itô Calculation

In this section we use the Itô calculus to obtain results on the overlap. We refer to [START_REF] Carmona | Strong disorder implies strong localization for directed polymers in a random environment ALEA Lat[END_REF] for the use of stochastic calculus in the study of the parabolic Anderson model.

Recall Z κ,β,t = E κ [exp{βH t (X)}]
and that

µ κ,β,t (f ) = E κ [f exp{βH t (X)}] Z κ,β,t . Note that d ln Z κ,β,t = 1 Z κ,β,t dZ κ,β,t - 1 2Z 2 κ,β,t d Z κ,β,t .
where

dZ κ,β,t =dE κ [exp{βH t (X)}] =E κ [d exp{βH t (X)}] =E κ β exp{βH t (X)} dW X(t) (t) + β 2 2 exp{βH t (X)} dt (36)
and therefore,

d Z κ,β,t = β 2 E ⊗2 κ 1 X(t)= X (t) exp{βH t (X)} exp{βH t ( X)} dt
with X an independent copy of X. Thus,

d ln Z κ,β,t = βµ κ,β,t (dW X(t) (t)) + β 2 2 1 -µ ⊗2 κ,β,t (X(t) = X(t)) dt
and upon integration we get

1 t ln Z κ,β,t = 1 t M t + β 2 2 1 - 1 t t 0 µ ⊗2 κ,β,s (X(s) = X(s))ds , ( 37
)
where M t is a square-integrable martingale with quadratic variation M given by

d M t dt = β 2 x∈Z d µ κ,β,t (X(t) = x) 2 = β 2 µ ⊗2 κ,β,t (X(t) = X(t)).
As a consequence we derive from (37) both the existence of the limit ( 14) and the relation [START_REF] König | A two cities theorem for the parabolic Anderson model[END_REF],

Ψ(κ, β) = β 2 2 1 -lim t→∞ 1 t t 0 µ ⊗2 κ,β,s X(s) = X(s)) ds , a.s.
Indeed,

d M t dt =β 2 x∈Z d µ κ,β,t (X(t) = x) 2 ≤β 2 sup y∈Z d µ κ,β,t (X(t) = y) x∈Z d µ κ,β,t (X(t) = x) ≤β 2 . ( 38
)
Thus, lim t→∞ 1 t M t = 0 and the rest is clear. Proof of Proposition 2.3. We have just proved statements (i). By (33), and since the complement of D κ is at most countable, we obtain the claims (ii).

Proof of Corollary 2.1. (i) is a consequence of ( 15) and of the definition (13) of Υ c .

(ii) follows from the integral formula (34). For instance, for 0

< β < β 1 such that (β 2 /2) -Ψ(κ, β) = (β 2 1 /2) -Ψ(κ, β 1 )
, the first term is smaller than the second one by monotonicity, and the difference

0 < (β 2 1 /2) -Ψ(κ, β 1 ) -(β 2 /2) + Ψ(κ, β) = β 1 β r 1 -J κ,r,∞ dr,
so that there exists r ∈ (β, β 1 ) such that J κ,r,∞ > 0. From this and similar considerations, we easily obtain the second statement. Under the assumption of the last claim, the equality (17) holds, proving the claim.

Proof of Proposition 2.2. The claim (i) is a consequence of ( 26) and [START_REF] König | A two cities theorem for the parabolic Anderson model[END_REF].The claim in (ii) is simply the second part of Proposition 3.1.

Proof of Theorem 2. Letting a T = 1 -I κ,β,T for a short notation, we have

a T := 1 - 1 T T 0 µ ⊗2 κ,β,t (X(t) = X(t))dt = 1 T T 0 x∈Z d µ κ,β,t (X(t) = x)µ κ,β,t (X(t) = x)dt = 1 T T 0 1 - x∈Z d µ κ,β,t (X(t) = x) 2 dt.
Let us denote

m t = µ κ,β,t (X(t) = x * (t)), b T = 1 T T 0 m t dt.
By splitting off the term for x = x * (t) from the sum of the terms for x = x * (t) in (39) we get,

a T = 1 T T 0   1 -(1 -m t ) 2 - x =x * (t) µ κ,β,t (X(t) = x) 2   dt = 2b T - 1 T T 0 m 2 t dt - 1 T T 0 x =x * (t) µ κ,β,t (X(t) = x) 2 dt.
Clearly, it implies that a T ≤ 2b T , but also an estimate in the reverse direction. Since

x =x * (t) µ κ,β,t (X(t) = x) 2 ≤ x =x * (t) µ κ,β,t (X(t) = x * (t)) × µ κ,β,t (X(t) = x) = (1 -m t )m t , by definition of x * , we obtain b T ≤ a T ≤ 2b T .
Therefore, (i) in Theorem 2 follows from points (i) of Corollary 2.1 and of Proposition2.2. We now turn to the proof of (ii). Repeating the same steps with µ κ,β,T instead of µ κ,β,t , except for the splitting according to x = y * T (t) or not, we see that

āT = 1 -E 1 T T 0 µ ⊗2 κ,β,T (X(t) = X(t))dt,
and bT = E 1 T T 0 µ κ,β,T (X(t) = y * T (t))dt,
are such that bT ≤ āT ≤ 2 bT .

Since āT = 1 -EJ κ,β,T , point (ii) in Theorem 2 follows from points (ii) of Corollary 2.1 and of Proposition 2.2.

Jump Distribution

Let E S denote the expectation for the simple, discrete time random walk S = {S(i); i ∈ N} on Z d with discrete time, and, for n ∈ N, T n,T = {(t 1 , . . . , t n ) : 0 < t 1 < . . . t n < T }. The quantity

E κ e βH T (X) |N(T, X) = n (39)
does not depend on κ and is equal to

(κT ) n n! e -κT -1 . . . T n,T E S n i=0 e β[W S(i) (t i+1 )-W S(i) (t i )] κe -κ(t i+1 -t i ) dt 1 . . . dt n = T -n n! . . . T n,T E S e β n i=0 [W S(i) (t i+1 )-W S(i) (t i )] dt 1 . . . dt n ,
where we have set t 0 = 0, t n+1 = T . Under the law on path space defined by (39), the jump times and jump values are independent, with respective distributions, uniform on [0, T ], and P S .

Proposition 5.1. The following limits exist a.s. and in L p , p ∈ [1, ∞), and are equal:

Λ(κ, β, r) = lim T →∞ T -1 ln E κ e βH T (X) ; N(T, X) = [rT ] (40)

= lim

T →∞,n/T →r T -1 ln E κ e βH T (X) ; N(T, X) = n

The limit is deterministic, convex in β, and concave in r.

We will use the following observation, which has been found to be useful in similar situations, where the subadditive (or superadditive) ergodic theorem does not apply. Lemma 5.1 (Stochastic superadditive lemma). Let U t be an integrable random process indexed by t in N or R + , such that (i) EU t+s ≥ EU t + EU s , s, t ≥ 0, and The claims now follows from (ii).

(ii) 1 t (U t -EU t ) * -→ 0 as t → ∞,
Proof of Proposition 5.1. For r ∈ Q + , we now check that Lemma 5.1 applies to the sequence U t = ln E κ e βHt(X) ; N(t, X) = rt , t ∈ N with rt ∈ N. First, define a probability measure on the set of paths D ∞ by µ

(r) κ,β,t (•) = µ κ,β,t (• |N(t, X) = rt), we have for s ∈ r -1 N, U t+s ≥ ln E κ e βH t+s (X) ; N(t + s, X) = r(t + s), N(t, X) = rt = ln E κ e βHt(X) e β[H t+s (X)-Ht(X)] ; N(t + s, X) = r(t + s), N(t, X) = rt Markov = U t + ln µ (r) κ,β,t exp U s • θ t,X(t) (θ time-space shift) Jensen ≥ U t + µ (r) κ,β,t U s • θ t,X(t) (41)
Since {W x : x ∈ Z d } are independent Brownian motions,

E U t + µ (r) κ,β,t U s • θ t,X(t) = EU t + Eµ (r) κ,β,t E U s • θ t,X(t) = EU t + EU s ,
which, together with (41), proves (i).

To show (ii) with a.s. convergence, we combine concentration and martingale inequalities. Let 1/2 < a < 1, and {T n } n≥1 be the sequence defined by T 1 = 1, T n+1 = T n + T a n . Then T n = n 1 1-a +o(1) as n → ∞. By Lemma 3.2 with A = {N(t, X) = rt}, it is easily seen that (ii) holds along the sequence T n using the Borel-Cantelli lemma. We now bridge the gaps. By Itô's formula,

U T = M T -(1/2) M T + β 2 T /2
for some continuous martingale M with (d/dt) M t ≤ β 2 for all t ≥ 0. Fix a sequence ǫ n → 0 with ǫ -1 n = n o (1) . For n large, β 2 T a n < ǫ n T n+1 , and then

Q sup Tn≤T ≤T n+1 |U T -U Tn -EU T +EU Tn | > 2ǫ n T n+1 ≤ E sup Tn≤T ≤T n+1 |M T -M Tn | > ǫ n T n+1 Doob ≤ (ǫ n T n+1 ) -2 E M T n+1 -M Tn ≤ (ǫ n T n+1 ) -2 β 2 (T n+1 -T n ),
which defines a summable series if we choose a ∈ (1/2, 1) large enough. By Borel-Cantelli, this completes the proof of (ii).

The limit Λ is convex in β as a limit of convex functions. We now check concavity in r. First note that V T,n = ln E κ e βH T (X) ; N(T, X) = n , which is equal to U T with r = T /n, satisfies

V T +T ′ ,n+n ′ ≥ ln E κ e βH T +T ′ (X) ; N(T, X) = n, N(T + T ′ , X) = n + n ′ = V T,n + ln µ (r) κ,β,T [exp V T +T ′ ,n+n ′ • θ T,X(T ) ] (r = T /n). (42)
Proceeding as in (41), and letting T, T ′ , n, n ′ → ∞ in such a way that n/T → r, n ′ /T ′ → r ′ , T /T ′ → λ/(1λ) with r, r ′ ≥ 0 and λ ∈ [0, 1], we get

Λ(κ, β, λr + (1 -λ)r ′ ≥ λΛ(κ, β, r) + (1 -λ)Λ κ, β, r ′ ),
i.e, the desired concavity.

Proof of Proposition 2.1. Writing the conditional expectation as a ratio, we now see that the limit lim

T →∞ T -1 ln E κ e βH T (X) |N(T, X) = [rT ]
exists, and is equal to Λ(κ, β, r)+I κ (r) by Cramér's theorem. This shows the existence of Γ, and also that Λ(κ, β, r) = Γ(β, r) -I κ (r) .

We now turn to the scaling relation. Under P κ (• |N(T, X) = n), X [0,T ] := (X t ; t ∈ [0, T ]) has n jumps on [0, T ] the values and times of which are independent, uniformly distributed. Then, with X (a) : s → X(s/a), the following image laws are equal

X (a) [0,aT ] • P κ x (• |N(T, X) = n) = X [0,aT ] • P κ x (• |N(aT, X) = n). (44)
Also, W

x (s) = a 1/2 W x (s/a) defines a collection, W (a) = {W (a) x : x ∈ Z d }, of independent standard Brownian motions. Denoting by t i , i = 1, . . . n, the jump times of X and t 0 = 0, t n+1 = T , we have by definition,

H W T (X) = n i=0 [W Xt i (t i+1 ) -W Xt i (t i )] = a -1/2 H W (a)
aT (X (a) ), and also,

1 T ln E κ e βH T (X) |N(T, X) = n = a aT ln E κ e βa -1/2 H W (a) aT (X (a) ) |N(T, X) = n = a aT ln E κ e βa -1/2 H W (a) aT (X) |N(aT, X) = n (45)
by (44). The first scaling relation follows from taking the limit T → ∞, n/T → r. The convexity in β follows from Hölder's inequality. For continuity we need to establish lim a→1 Γ(β, ar) = Γ(β, r). But, by scaling, Γ(β, ar) = aΓ(β/ √ a, r) and the result follows from continuity in β.

Here is a direct consequence of (43) and of the scaling.

Corollary 5.1. For all r, Λ is jointly convex in (κ, β), and Λ(κ, β, r) = aΛ(a -1 κ, a -1/2 β, a -1 r).

(46)

Proof of Proposition 2.2. For a, b > 0, write

E κ e βH T (X) =
aT ≤n≤bT E κ e βH T (X) ; N(T, X) = n + E κ e βH T (X) ; N(T, X) > bT +E κ e βH T (X) ; N(T, X) < aT .

(

) 47
By the concentration inequality in Lemma 3.2, as T → ∞, X) ; N(T, X) > aT -E ln E κ e βH T (X) ; N(T, X) > aT → 0 a.s. and in L p . By Jensen's inequality,

T -1 ln E κ e βH T (
T -1 E ln E κ e βH T (X) ; N(T, X) > bT ≤ β 2 /2 -I κ (B),
which can be made arbitrarily negative by taking b large. Similarly,

T -1 E ln E κ e βH T (X) ; N(T, X) < aT ≤ β 2 /2 -I κ (a),
which can be made arbitrarily negative by taking a small. Thus, for a sufficiently small and b sufficiently large, lim

T →∞ T -1 ln E κ e βH T (X) = lim T →∞
T -1 ln E κ e βH T (X) ; aT < N(T, X) < bT .

Define Γ T (β, r) = T -1 E ln E κ e βH T (X) |N(X, T ) = [rT ] which is a convex function of β converging point-wise to Γ(β, r) which is also convex in β. The conditional version of the concentration inequality holds due to cancellation, that is

P 1 T ln E κ e βH T (X) |N(T, X) = [rT ] -Γ T (β, r) ≥ u ≤ 2 exp - u 2 T 2β 2 . Write r = sa with s ranging over C ≡ [1, b a ] ∩ n -1 Z.
Since the number of points in C grows like n, we conclude by Borel-Cantelli that

P sup s∈C 1 n ln E κ e βHn(X) |N(n, X) = [saT ] -Γ n (β, sA) > ǫ/4 i.o. = 0.
We also note that by convexity of both Γ T and Γ that Γ T (β, r) converges uniformly for β in a compact interval to Γ(β, r). By (45) it follows that Γ T (β, ar) = aΓ aT (a -1/2 β, r) which implies uniform convergence of Γ T (β, r) to Γ(β, r) for r in a compact interval.

Finally, for 0 < a < b < ∞, r∈[a,b]∩n -1 Z E κ e βHn(X) ; N(n, X) = [rn] = r∈[a,b]∩n -1 Z E κ [exp{βH n (X)}|N(n, X) = [rn]] ×P κ (N(n, X) = [rn]) = e o(n) r∈[a,b]∩n -1 Z e Γn(β,r)n P κ (N(n, X) = [rn]) = e o(n) r∈[a,b]∩n -1 Z e Γ(β,r)n P κ (N(n, X) = [rn]) = e o(n) r∈[a,b]∩n -1 Z e (Γ(β,r)-Iκ(r))n . Thus, T -1 ln E κ e βH T (X) ∼ T -1 ln E κ e βH T (X) ; aT < N(T, X) < bT = o(1) + T -1 ln r∈[a,b]∩T -1 Z e (Γ(β,r)-Iκ(r))T = sup Γ(β, r) -I κ (r); r ∈ [a, b] + o(1) → Ψ(κ, β), T → ∞, (48)
by the standard Laplace method.

Remark 5.1. We have shown that Ψ(κ, β) = sup{Λ(κ, β, r); r ≥ 0}.

Since I κ is non negative and zero if r = κ, and by (43), we have Γ(β, r) = sup{Λ(κ, β, r); κ > 0}.

Proof of Theorem 1. The first statement directly follows from the previous results. The large deviation principle ( 12) is proved in a manner similar to Proposition 2.2

Proof of Proposition 2.5. By [START_REF] Carmona | Strong disorder implies strong localization for directed polymers in a random environment ALEA Lat[END_REF],

Γ(β, r) β = r Γ(r -1/2 β, 1) β = √ r Γ(r -1/2 β, 1) r -1/2 β . (49)
We claim lim

β→∞ β -1 Γ(β, 1) = α. (50)
i.e., [START_REF] Seppäläinen | Bounds for scaling exponents for a 1 + 1 dimensional directed polymer in a Brownian environment[END_REF] holds for r = 1, and that for all β, β -1 Γ(β, 1) ≤ α. Indeed,

1 β Γ(β, 1) = lim T →∞ 1 βT ln E κ [e βH T (x) |N(T, X) = T ] ≤ lim T →∞ 1 βT βA T =α. (51)
Following the developments in [START_REF] Cranston | Lyapunov exponents for the parabolic Anderson model[END_REF], given ǫ > 0, one can find a path γ ∈ D T with N(T, γ) = T and the jump times of γ are separated by ǫ and H T (γ) > (αδ)T where δ → 0 as ǫ → 0. Moreover, writing η ∼ γ to mean that the jump times of η are within ǫ/3 of the jump times of γ and the two paths jump to the same sites at these jump times we have

P (X ∼ γ|N(T, X) = T ) ≥ T ! T T ǫ 6 T ∼ √ T ǫ 6e
T .

In addition, for X ∼ γ one has eventually, H T (X) ≥ (αδ)T with δ → 0 as ǫ → 0. Thus,

lim β→∞ lim T →∞ 1 βT ln E κ [e βH T (x) |N(T, X) = T ] ≥ lim β→∞ lim T →∞ 1 βT ln E κ [e βH T (x) , X ∼ γ|N(T, X) = T ] ≥ lim β→∞ lim T →∞ 1 βT ln e β(α-δ)T P κ (X ∼ γ|N(T, X) = T ) ≥ lim β→∞ lim T →∞ 1 βT ln e β(α-δ)T √ T ǫ 6e T ≥(α -δ) (52)
and letting ǫ → 0 and therefore δ → 0 we have established the claim (50). Then, (50) and ( 52) imply (i), uniformly for r ∈ [ǫ, ǫ -1 ] for all ǫ ∈ (0, 1]. By ( 49) and (51) we have Γ

(β, r) ≤ αβ √ r, (53)
Combining ( 9) and (53),

Ψ(κ, β) ≤ sup r≥0 αβ √ r -I κ (r) .
For the converse direction, consider r m the (unique) maximizer of αβ [START_REF] Rovira | On the Brownian-directed polymer in a Gaussian random environment[END_REF].

√ r -I κ (r). Note that r m ∼ α 2 β 2 4 ln 2 (β 2 /κ) as β 2 /κ → ∞ and write Ψ(κ, β) ≥ Γ(β, r m ) -I κ (r m ) = r m Γ(r -1/2 m β, 1) -I κ (r m ) (by (7)) = r m -1/2 β(1 + ǫ(r -1/2 m β)) -I κ (r m ) (lim u→0 ǫ(u) = 0, see (50)) = (1 + ǫ(r -1/2 m β)) × r.h.s. of
Thus, the first claim in (ii) is proved. The second one is clear. We finish by (iii), which means that all maximizers are of the indicated order of magnitude. For ǫ > 0, using for some positive δ, δ ′ .

Regularity of the favourite attributes

In this section, we give the Proof of Proposition 2.4. We denote by x = max i≤d |x i | the supremum norm on the lattice. Let 0 < a < b < ln 2.

(i) Observe that provided n is sufficiently large that 2 n /t > κ, P sup

2 n ≤ x <2 n+1 f (x,
2 n ≤ x <2 n+1 f (x, t) > e -bn2 n ) < ∞,
and consequently, given almost any realization of W, there is an N 1 such that f (x, t) ≤ e -bn2 n , for x ≥ 2 n , n ≥ N 1 .

This implies that for some N 1 = N 1 (W), (55)

In the same way as before, we now conclude that f (x, t, T ) exhibits superexponential decay in x and as well that for some N 2 = N 2 (W, T ), sup

x∈Z d f (x, t, T ) = max x∈Z d , x ≤2 N 2
f (x, t, T ), 0 ≤ t ≤ T.

(iii) For X ∈ D t , define the time-reversed path X(t) ∈ D t with X(t) (s) = X(ts). Since the symmetric simple random walk is reversible, we have f (x, t) = E κ [exp{βH t (X)}δ x (X(t))] = E x κ [exp{βH t ( X(t) )}δ 0 (X(t))]. By [START_REF] Amir | Quastel: Probability Distribution of the Free Energy of the Continuum Directed Random Polymer in 1+1 dimensions[END_REF], we obtain f (x, t) = δ x (0) + κ t 0 ∆f (x, s)ds + β t 0 f (x, s) • dW x (s).

Since f is positive, we see that f (x, •) has the same regularity as W x .

(iv) Define    H x t,T (X) = y =x T t δ y (X(st))[dW y (s) -dW x (s)], H x t (X) = H x 0,t (X), f x (x, t, T ) = f (z, t, T ) × e -βWx(T ) .

(56) By Markov property, we have f x (z, t, T ) = f x (z, t) × g x (z, t, T ), (57) with f x (z, t) = E κ [e βH x t (X) δ z (X(t)), g x (z, t, T ) = E z κ [e βH x t,T (X) ].

By reversibility under P κ , f x (z, t) = E z κ [e βH x t ( X(t) ) δ 0 (X(t))].

Similar to (1), we obtain for z = x,

By definition of g x (z, t, T ), we have a similar identity, g x (x, t, T ) = 1 + κ T t ∆ z g x (x, s, T )ds.

Combining this with (56), ( 57), (58), we conclude that t → f x (z, t, T ) is differentiable on (0, T ) with derivative d dt f (x, t, T ) = e βWx(T ) κ g x (x, t, T )∆ z f x (x, t)f x (x, t)∆ z g x (x, t, T ) , which is continuous this interval, with probability 1.

Proposition 2 . 5 .

 25 (i) Almost surely, lim β→∞ β -1 Γ(β, r) = α √ rlocally uniformly for r ∈ (0, +∞).(ii) As β 2 /κ → ∞,

  53) and the definition ofr m , sup{Γ(β, r) -I κ (r) : |rr m | ≥ r m ǫ} ≤ sup αβ √ r -I κ (r); |rr m | ≥ r m ǫ ≤ (1δ) sup{αβ √ r -I κ (r); r ≥ 0} ≤ (1δ ′ )Ψ(κ, β)

  sup x∈Z d f (x, t) = max x∈Z d , x ≤2 N 1 f (x, t),which implies the desired property.(ii) Setting τ x = inf{s ≥ 0 : X(s) = x}, we have P ( sup2 n ≤ x <4 d 2 n+1 sup t≤T f (x, t, T ) > e -bn2 n ) ≤ 4 d 2 nd max 2 n ≤ x <2 n+1 P (sup t≤T f (x, t, T ) > e -bn2 n ) ≤ 4 d 2 nd e bn2 n max 2 n ≤ x <2 n+1 E sup t≤T E κ [exp{βH T (X)}δ x (X(t))] ≤ 4 d 2 nd e bn2 n max 2 n ≤ x <2 n+1 EE κ sup t≤T exp{βH T (X)}δ x (X(t)) = 4 d 2 n e bn2 n max 2 n ≤ x <2 n+1 EE κ [exp{βH T (X)}; τ x ≤ T ] = 4 d 2 nd e bn2 n max 2 n ≤ x <2 n+1 E κ [E [exp{βH T (X)}]] ; τ x ≤ T ] ≤ 4 d 2 nd e bn2 n e β 2 T /2 P κ (N(T, X) ≥ 2 n ) ≤ 4 d2 nd e bn2 n e β 2 T /2 e -2 n log(2 n /κT )+(2 n /T -κ)T .

0 ∆ 0 ∆

 00 f x (z, t) = δ z (0) + κ t z f x (z, s)ds + β t 0 f x (z, s) • d[W z (s) -W x (s)],though the last term vanishes for z = x, f x (x, t) = δ x (0) + κ t z f x (x, s)ds.

  where * -→ is some stochastic mode of convergence (a.s., in probability, in L p , . . . ).

	Proof of Lemma 5.1. From (i) and the superadditive lemma,
		lim t→∞	EU t t	= sup	EU t t	; t ≥ 0 .
	Then, as t → ∞,	U t t	* -→ sup	EU t t	; t ≥ 0
	with the same mode of convergence.

  t) > e -bn2 n ≤ 4 d 2 nd max 2 n ≤ x <2 n+1 P (f (x, t) > e -bn2 n ) ≤ 4 d 2 nd e bn2 n max 2 n ≤ x <2 n+1 EE κ [exp{βH t (X)}δ x (X(t))] = 4 d 2 nd e bn2 n max 2 n ≤ x <2 n+1 E κ [E [exp{βH t (X)}]] δ x (X(t))] ≤ 4 d 2 nd e bn2 n e β 2 t/2 P κ (N(t, X) ≥ 2 n ) ≤ 4 d 2 nd e bn2 n e β 2 t/2 e -2 n log(2 n /κt)+(2 n /t-κ)t .

			(54)
	Thus,	∞
		P (	sup
		n=0

Acknowledgements: FC thanks the Fields institute for hospitality during this work.

Research of second author supported by a grant from NSF, DMS 0854940.


Karim Belharet

Member, IEEE David Folio

Member, IEEE Antoine Ferreira

Simulation and Planning of a Magnetically Actuated Microrobot Navigating in Arteries

Keywords: Microrobotics, minimally invasive surgery, blood flow simulation, anisotropic path planning

published or not. The documents may come

I. INTRODUCTION

P REOPERATIVE surgical simulation and planning is a critical stage in medical decision-making and in the design of novel microrobotic minimally invasive surgery (MIS) techniques [START_REF] Nelson | Microrobots for minimally invasive medicine[END_REF][START_REF] Miloro | Removing vascular obstructions: a challenge, yet an opportunity for interventional microdevices[END_REF]. Due to the lack of anatomic consistency, preoperative planning with virtual and physical models of the system could improve the chances of achieving predictable intraoperative results. To this aim, it is mandatory to cope with potential conflicts between the intervention purposes of the microrobotic surgical system and the biological laws governing the patient body. To improve microrobotic MIS procedure significant progress has been made in three dimensional (3D) imagers [START_REF] Satava | Virtual reality, telesurgery, and the new world order of medicine[END_REF]; modeling software for anatomy and vascular networks representation [START_REF] Sun | Local morphology fitting active contour for automatic vascular segmentation[END_REF]; and biomechanical modeling of interactions between tissues and MIS tools. These advances have been mainly developed for catheter or guidewire-based operations [START_REF] Ganji | Catheter kinematics for intracardiac navigation[END_REF][START_REF] Tang | A stable and real-time nonlinear elastic approach to simulating guide-wire and catheter insertions based on cosserat rod[END_REF], or gastrointestinal endoscopic capsules [START_REF] Carpi | Magnetically controllable gastrointestinal steering of video capsules[END_REF][START_REF] Ciuti | A comparative evaluation of control interfaces for a robotic-aided endoscopic capsule platform[END_REF]. In contrast, there has been fewer works in the area of untethered microrobots for cardiovascular interventions [START_REF] Nelson | Microrobots for minimally invasive medicine[END_REF][START_REF] Tamaz | Real-time MRI-based control of a ferromagnetic core for endovascular navigation[END_REF]. To date, even in the case of a successfully performed procedure, the choice of a microrobot's design (shape, size, materials, locomotion abilities. . . ) plays a key role. Indeed, these factors affect the microrobot propulsion that depends mainly on the motile, fluidic and friction forces between the endovascular device and the vessel. Therefore, the development of training applications or interactive planning systems, where the robotician or physician is able to design different microrobots and then test their behaviors in a patient-specific context, will be very helpful. A few applications focusing on interventional MIS have been developed or commercialized [10,[START_REF]Sensei(r) x robotic catheter system[END_REF] until now. These applications include interactive models of catheters [START_REF] Ganji | Catheter kinematics for intracardiac navigation[END_REF], endoscopic capsules [START_REF] Carpi | Magnetically controllable gastrointestinal steering of video capsules[END_REF], but do not address the challenging topic of modeling untethered microrobots [START_REF] Pawashe | Modeling and experimental characterization of an untethered magnetic micro-robot[END_REF] and their complex interactions with the blood flow and vessels walls [START_REF] Nelson | Microrobots for minimally invasive medicine[END_REF]. Most proposed solutions use finite element modeling (FEM)-based representations, or computational fluid dynamics (CFD) simulations. Such approaches are suitable when accurate modeling or simulation are a primary concern. However, they usually rely on commercial software (such as Fluents or Comsol), and the computation time (usually in hours) is incompatible with interactive simulation or clinical practice. Moreover, planning a reliable and feasible path against the blood flow and vessel networks constraints is an important issue to be addressed for innovative microrobotaided MIS applications. For instance, S. Martel et al. propose to use potential fields and the breadth-first search algorithm to find a pre-planned path [START_REF] Sabra | Use of 3D potential field and an enhanced breadth-first search algorithms for the path planning of microdevices propelled in the cardiovascular system[END_REF] which are then delimited by waypoints [START_REF] Tamaz | Real-time MRI-based control of a ferromagnetic core for endovascular navigation[END_REF]. A first drawback in such discrete graphsearch algorithms is that they suffer from metrication errors [START_REF] Pêtrès | Path planning for autonomous underwater vehicles[END_REF]. Moreover, it is difficult to cope explicitly with the anisotropy. Finally, most untethered microrobot propulsion schemes based on magnetic pulling have to face important constraints related to coils technology. The planned path has to reduce the maximum pulling gradients being applied to the coils, the operation time and to maintain the rising temperature of the coils within operating limits. It actually impacts different endovascular magnetic actuation systems based on Magnetic Resonance Navigation (MRN) [START_REF] Tamaz | Real-time MRI-based control of a ferromagnetic core for endovascular navigation[END_REF], or Octomag system [START_REF] Kummer | Octomag: An electromagnetic system for 5-dof wireless micromanipulation[END_REF]. Thus, the planned pathway must satisfy these considerations to be reliable.

This work aims to propose a preoperative microrobotic surgical planning application with blood flow simulation that may use patient-specific data. In particular, the computeraided optimization of the untethered microrobots for innovative MIS procedure is the main focus of our work. First, in Section II, after introducing the system model used by the microrobotic surgical simulation, we present the surgical planning framework. The proposed framework is based on the Fast Marching Method (FMM) [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF]. In our previous work [START_REF] Belharet | Endovascular navigation of a ferromagnetic microrobot using MRI-based predictive control[END_REF][START_REF]3D MRI-based predictive control of a ferromagnetic microrobot navigating in blood vessels[END_REF], the classical FMM with isotropic cost function has been proposed to extract the vessel centerline path. This first solution is suitable mainly to navigate far away from the vessel wall. However, as navigation in real vascular networks is directionally constrained by the biological environment, a new framework based on the anisotropic FMM [START_REF] Sethian | Ordered upwind methods for static hamilton-jacobi equations: theory and algorithms[END_REF] is proposed in this paper. Furthermore, the proposed method allows a multi-parameter optimization of the navigation path that takes into account: i) the technology of the magnetic system; ii) the microrobot design; and iii) the physiology of the patient. Finally, our approach ensures a smooth conveyance of the microrobot to destination robust to cardiovascular navigation constraints while minimizing the energy expenditure, to reduce the gradient pulling and system overheating. Specifically, we aim to get a feasible path by finding the min-energy pathway (which is not necessarily the shortest). Section III briefly describes our developed microrobotic surgical simulation and planning application, before illustrating its performance. As an illustration, we consider the treatment of peripheral arterial diseases (PAD), such as artery narrowing (stenosis). For instance, one possible surgical application is to target atherosclerotic plaques in stenosed arteries. One of the most common strategies is angioplasty, occasionally with stent implantation. Despite the high success rate, these solutions are not satisfactory, as vessel damage, infection and restenosis can appear at a later stage. Otherwise, a surgical microrobot for such vascular applications offers some novel solutions [START_REF] Nelson | Microrobots for minimally invasive medicine[END_REF]. Specifically, microrobotic intervention based on chemical, mechanical or mechano-chemical could be achieved to treat vascular obstructions [START_REF] Miloro | Removing vascular obstructions: a challenge, yet an opportunity for interventional microdevices[END_REF]. Finally, we discuss some open issues, evaluate the computational efficiency, and outline some future extends in Section IV. This study is concluded in Section V.

II. METHODS

A. A modeling overview

Different types of devices could be used for MIS interventions from catheter and guidewires, to wireless microcapsules or microcarriers [START_REF] Nelson | Microrobots for minimally invasive medicine[END_REF]. In this framework, any kind of untethered microdevice (termed as microrobot throughout the text) that has to perform endovascular navigation within blood flow is considered. More precisely, we consider any microsystem that could be modeled as:

m v = f a (a) + f d (v) + f Ext (x) (1)
where v is the velocity of the device; m its mass; f a (a) is the controlled force, that is related to the control action a ∈ A, with A the admissible controls set; f d (v) denotes the fluid flow hydrodynamic drag force; and f Ext (x) are all other external forces, as depicted Fig. 1. Different external forces could be added to the model (1), such as the weight (f g ), electrostatic (f e ), van der Waals (f v ), contact or steric forces. The interested reader may refer to [START_REF] Arcèse | Endovascular Magnetically-Guided Robots: Navigation Modeling and Optimization[END_REF] for detailed formulation of microforces. Fig. 2 illustrates the evolution of these microforces with varying microrobot radius r and wrt. its distance to vessel wall δ. In particular, Fig. 2(f) shows that close to the centerline of the vessel (δ = R) the electrostatic and the van der Waals microforces are negligible compared to the other relevant forces. In addition, if the microrobot is never in contact with the vessel wall, the contact force could be omitted. Finally, most of the time because of the size of the microsphere, the effect of Brownian motion is also neglected.

B. Navigation Planning in Flow

1) Minimal Path Planning: Classically, a navigation pathway could be defined as a curve P linking a starting point x 0 to any point x. In the general case, a curve is a minimal path, also called a geodesic, wrt. the metric ϖ, if it globally minimizes the energy functional:

E ϖ (P) = ∫ P ϖ(P(l), Ṗ(l))dl, (2)
with l the arclength parameter. The solution to the minimum path finding problem (2) could be obtained through the computation of the distance map U : C → R + , defined as: U(x) = min P E(P), for any P linking x to x 0 in the domain C ⊂ R d (with d = 2, 3, 4. . . ). In this context, Sethian [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF] has proposed the Fast Marching Method (FMM). Indeed, the FMM converges to a smooth solution in the continuous domain even when it is implemented on a sampled environment, contrary to discrete graph-search algorithms. Especially, in the isotropic case, the FMM satisfies the Eikonal equation:

{ ∥∇U(x)∥ = w(x), ∀x ∈ C U(x 0 ) = 0 (3)
where the isotropic cost function w : x ∈ C → w(x) ∈ R + is the metric. As the map U has a single local minimum, the geodesic can be retrieved with a simple gradient descent on U from a targeted seed x g to x 0 . The FMM provides a continuous solution to the minimum path problem by employing upwind differences and a causality condition [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF]. The key issue of the FMM is then to get an appropriate metric w(x) which drives the front expansion efficiently to find a geodesic P. In previous works [START_REF] Belharet | Endovascular navigation of a ferromagnetic microrobot using MRI-based predictive control[END_REF][START_REF]3D MRI-based predictive control of a ferromagnetic microrobot navigating in blood vessels[END_REF], we have proposed to design such a metric on spatial consideration:

w : x ∈ C → Vesselness(x) ∈ R + (4)
where Vesselness(x) is a vesselness enhancement function.

A typical Vesselness(x) function could be designed using a multi-scale analysis of the Hessian matrix, as proposed by Frangi [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] or Sato [START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF]. This approach then allows to find the vessel centerline geodesic P c , as the considered vesselness filter gives their maximal response in the vessel center. We have applied this procedure to different sets of data, in 2D [START_REF] Belharet | Endovascular navigation of a ferromagnetic microrobot using MRI-based predictive control[END_REF] as well as in 3D [START_REF]3D MRI-based predictive control of a ferromagnetic microrobot navigating in blood vessels[END_REF].

As mentioned previously the contact, electrostatic and van der Waals microforces are negligible when the microrobot navigates close to the centerline P c (cf. Fig. 2). P c is then the best solution to navigate far away from the vessels wall. However, navigation in real vascular networks is constrained by the biological environment, and can strongly modify the system behavior, leading to incompleteness in real intraoperative applications.

2) Anisotropic Path Planning: Navigating in blood flow implies some directional constraints. When the path finding problem has some preferred directions of travel it becomes anisotropic, i.e. the solutions depend on the configuration and directional constraints. In [START_REF] Sethian | Ordered upwind methods for static hamilton-jacobi equations: theory and algorithms[END_REF], the authors noticed that the FMM could be applied to isotropic or anisotropic problems, and they provided a numerical scheme to solve static Hamilton-Jacobi equations of the form:

{ ∥∇U(x)∥ = ϖ(x, ∇U (x) ∥∇U (x)∥ ), ∀x ∈ C U(x 0 ) = 0 (5)
where the anisotropy expansion is the result of the dependence of ϖ on ∇U(x). Nevertheless, the computational complexity of the FMM in the anisotropic case is growing significantly with the amount of anisotropy. The literature provides different variants of the FMM to deal with anisotropy, such as using a recursive approach [START_REF] Konukoglu | A recursive anisotropic fast marching approach to reaction diffusion equation: Application to tumor growth modeling[END_REF], or iterative schemes [START_REF] Bornemann | Finite-element discretization of static hamilton-jacobi equations based on a local variational principle[END_REF]. Furthermore, in [START_REF] Pêtrès | Path planning for autonomous underwater vehicles[END_REF] anisotropic cost functions have been also proposed. The proposed basic idea is to take advantage of the original FMM efficiency, and to model the anisotropy thanks to a vector f of a field of force F, leading to the following anisotropic metric:

ϖ : (x, f ) ∈ C × F → ϖ(x, f ) ∈ R + (6)
Such an anisotropic cost function leads to anisotropic frontpropagation, and gives better geodesic than isotropic FMM. To embed the directional constraints f within such an anisotropic FMM different specific cost functions have been proposed mainly for autonomous underwater vehicles (AUVs) [START_REF] Pêtrès | Path planning for autonomous underwater vehicles[END_REF][START_REF] Soulignac | Feasible and optimal path planning in strong current fields[END_REF]. The drag and flow forces have a significant effect on the motion of the AUV, which is a common problem with our microrobot planning problem. However, the proposed solution mainly deals with either the min-distance or the min-time path finding problem, and not minimizing directly neither the dynamics nor the motion applied to the system. Moreover, most existing approaches classically simplify the problem by assuming constant vehicle speed, or only spatially varying flow fields without time variation. Such limitations can generate physically unfeasible paths. Especially, in the presence of strong flows, incompleteness or incorrectness issues are reported [START_REF] Soulignac | Feasible and optimal path planning in strong current fields[END_REF].

3) Minimum-feasible-energy Planning: Most microrobot propulsion schemes usually require a limitation of the power consumption, and have a maximum available power. Moreover, considering only either the min-time or min-distance problem usually leads to compute high motion forces, that is reach the goal as fast as possible. But in our context, high motion implies high drag forces, and then high motion force f a (limited by the maximum available magnetic forces) to be able to reach the targeted area. To minimize the energy expenditure, a basic idea is to use a classical power functional expression, that is:

p(x, a) = ⟨f (x, a) • v(x)⟩ (7)
where ⟨•⟩ is the inner product of R d ; v is a velocity; and f (x, a) is the applied resultant force that embeds the controlled action a ∈ A. Using the optimal control approach, we propose the following anisotropic cost function:

ϖ(x, f ) = max a∈A ⟨f (x, a) • v(x)⟩ (8)
This anisotropic cost formulation, as in [START_REF] Pêtrès | Path planning for autonomous underwater vehicles[END_REF][START_REF] Soulignac | Feasible and optimal path planning in strong current fields[END_REF], naturally incites the device to move in the direction of the force f . Especially, ϖ(x, f ) will imply that the FMM propagation [START_REF] Ganji | Catheter kinematics for intracardiac navigation[END_REF] follows the force field constraints, that is ∇U (x) ∥∇U (x)∥ ∝ f . Finally, this new anisotropic formulation (8) using optimal control aspects embeds the system motion wrt. a finite defined admissible control set A. This ensures that the given geodesic P is really feasible by the system and minimizes energy expenditure. Indeed, the metric ϖ(x, f ) allows to consider both the biological constraints and the microrobot capabilities.

III. RESULTS

A. Simulation and Planning Tool

The proposed framework has been integrated in a dedicated interactive microrobot surgical simulation and planning application. This application, depicted on Fig. 3, is based on the Visualization Toolkit (VTK) [START_REF]The visualization toolkit (VTK)[END_REF], the Insight Segmentation and Registration Toolkit (ITK) [START_REF]Insight segmentation and registration toolkit (ITK)[END_REF], and the Qt Toolkit [START_REF]Qt-a cross-platform application and UI framework[END_REF]. Through the user interface (Fig. 3(b)), the surgeon realizes the following operating mode iteratively: a) Specifications: The surgeon provides different parameters in the software (in the left panel on Fig. 3(b)) related to: i) the physiology of the patient (blood flow and viscosity, vessel's size and topology); ii) the technology of the magnetic system (coil's slew rate, maximum magnetic gradients, power limitations); and iii) the microrobot design (size, density, magnetic saturation).

b) Simulations: The surgeon initiates the simulation and planning process wrt. the specification constraints.

c) Data analysis: The surgeon analyses the power functional mapping and force field wrt. specifications. This iterative process is conducted until the preoperative planning of plaque removal specifications are satisfied. We focus in the on the main contributions of this paper that is the simulation and reliable vascular navigation planning.

B. Case study: peripheral arterial diseases

To assess the overall framework, we have performed a series of tests for both the microrobot within blood vessel model simulation and preoperative microrobotic MIS planning. We illustrate the proposed application for the treatment of peripheral arterial diseases (PAD). One possible application is to locate atherosclerotic plaque in stenosed arteries. In particular, stenosis implies an abnormal narrowing of the vessel which causes a decrease in blood flow. Different plaque removal tasks have already been realized by magnetically actuated microrobots navigating in occluded blood vessels [START_REF] Miloro | Removing vascular obstructions: a challenge, yet an opportunity for interventional microdevices[END_REF]. In [START_REF] Korin | Shear-activated nanotherapeutics for drug targeting to obstructed blood vessels[END_REF], shear-targeting of a thrombolytic drug in a mouse arterial thrombosis model using a microcarrier has been successfully tested. In [START_REF] Yu | Novel electromagnetic actuation system for three-dimensional locomotion and drilling of intravascular microrobot[END_REF], the rotation of the magnetic microrobot is at the origin of a drilling action through an occlusion in a vessel. In the following, we propose a preoperative planning of a magnetic microrobot navigating within the vasculature to the stenosed artery for plaque removal. We will consider here as case study the real Magnetic Resonance Angiogram (MRA) clinical dataset depicted on Fig. 4. As the proposed framework works in C ⊂ R d (with d = 2, 3. . . ), for the sake of clarity of the representation, we will only show here a Maximum Intensity Projection (MIP) image of the MRA. Here, the patient has a single level disease represented by an isolated stenosis in the iliac artery. As it is still a challenging issue to identify automatically the patient pathology and its location, the PAD is manually defined by the user.

C. Simulation Results

As previously mentioned, the simulation is based on a physics-based model (presented in details in [START_REF] Arcèse | Endovascular Magnetically-Guided Robots: Navigation Modeling and Optimization[END_REF]). When interactive simulation is required (for training, intervention planning, etc.), a trade-off between computation time and accuracy has to be found. In our application, computational efficiency is more important than small-scale details, as we aim to increase the users interaction capabilities, whereas precision could be achieved using FEM or CFD simulation. To speed-up either the blood flow velocity or the force fields simulation computation, the basic FMM is extended to be able to compute these vector fields during the front propagation (see Fig. 3(a)). Indeed, using Vesselness(x) isotropic cost function, the FMM front will extend only within the vascular structure. This allows to take advantage of the FMM efficiency, without using the whole dataset grids, and significantly reduce the overall computational time.

The blood flow velocity field simulation of the considered case study is illustrated on Fig. 5. The interactive simulation tool takes into account that blood flow velocity intensity v f on the vessel centerline decreases wrt. vessels radii, that is from about ∥v ad f ∥ ≈ 400 mm/s in descending aorta with a radius R ≈ 6.5 mm, to ∥v ei f ∥ ≈ 250 mm/s in external iliac artery with a radius R ≈ 3.5 mm. As one can see, the blood flow velocity field is modified in the considered stenosis case (narrowed vessel on Fig. 5). The corresponding force field that will act on a neodymium microrobot with a radius r = 500 µm is shown on Fig. 6. In this simulation, only the weight f g (in the -→ g -axis direction), drag f d and electrostatic f e forces are calculated. Let us notice that any force that can have an analytical expression and can be embedded in [START_REF] Nelson | Microrobots for minimally invasive medicine[END_REF] can easily be added to the application. Finally, the realized interactive simulation application is also able to deal with periodic time-varying pulsatile blood flows, as presented on Fig. 7. The user defines the periodic pulsatile blood flow parameters [START_REF] Arcèse | Endovascular Magnetically-Guided Robots: Navigation Modeling and Optimization[END_REF], a number of time-steps, and can then browse the time-frame of either the blood flow velocity or the applied force field. Fig. 7(a) illustrates a typical periodic time-varying pulsatile blood flow that takes place in the descending aorta, and the row 7(b) describes the corresponding velocity field. Fig. 7(c) and 7(d) show the resulting force fields applied on the microrobot for two case studies, a microrobot of radius of r = 250 µm and r = 1.5 mm, respectively. This demonstrates that when the microrobot's radius r or weight f g increases, the resultant force becomes more important.

D. Microrobotic surgical planning demonstration

Once the microrobot is specified (size, materials, etc.), the system's planning step only requires as user's input the starting x 0 and ending points x g . Next, to deal with the anisotropy of the media, the blood flow and force fields simulation results are used to compute the power functional map [START_REF] Carpi | Magnetically controllable gastrointestinal steering of video capsules[END_REF]. The key issue is then to generate the controlled force field f a (a), ∀a ∈ A. First, the action magnitude has to be defined wrt. available actuation command. The action magnitude is usually the controlled magnetic gradient (a = ∇B), which implies the following magnetic force on a spherical device:

f a (∇B) = Dτ m 4 3 πr 3 (M.∇) ∇B, ( 9
)
where r is the microsphere radius; τ m ∈ [0; 1] is the volume ratio of the magnetic material; M its magnetization; and D ∈ [0; 1] is the duty-cycle, i.e., the ratio between the duration of the applied magnetic gradient for propulsion to the total period T . The magnetic gradient used to propel a magnetic microrobot is assumed to be generated from various medical magnetic steering systems: ElectroMagnetic based Actuation (EMA) systems such as the Magnetic Catheter Steering System (a larger version of the Octomag system [START_REF] Kummer | Octomag: An electromagnetic system for 5-dof wireless micromanipulation[END_REF]) or the one developed in [START_REF] Yu | Novel electromagnetic actuation system for three-dimensional locomotion and drilling of intravascular microrobot[END_REF]; and Magnetic Resonance Navigation (MRN). The EMA systems are capable of generating maximum magnetic gradient fields up to a max = 350 mT/m in closed-loop navigation control. In this case, the duty cycle is considered equal to D = 1 since the imaging modality is provided by either optical microscopy or X-ray imaging. These systems are more suitable for the navigation of millimeter-sized hard magnetic microrobots. The MRN system is based on clinical MRI systems, without any hardware upgrades. Such magnetic systems are usually limited to maximum values to a max = 80 mT/m in closed-loop navigation control. Generally, the duty cycle is D < 1, since MR-imaging sequences and MRtracking are provided by the same scanner to gather feedback positional information during MRN [START_REF] Tamaz | Real-time MRI-based control of a ferromagnetic core for endovascular navigation[END_REF]. The MRN approach is more suitable for micrometer-sized soft magnetic microrobots. However, it should be noticed that experimental results have shown that therapeutic particles (known as TMMC), 50-60 µm diameter were navigated in arteries (rabbit models) and that MRN was done with gradient coils (capable of 450 mT/m with inner diameter suitable for limbs or head interventions) installed in the tunnel of clinical MRI scanners [START_REF] Pouponneau | Coencapsulation of magnetic nanoparticles and doxorubicin into biodegradable microcarriers for deep tissue targeting by vascular MRI navigation[END_REF].

Our proposed FMM strategy can be adapted to each specific case by taking into account the duty-cycle D, the slew rate of the magnetic coils, and the maximum temperature coils. Secondly, the action direction should be defined. To this aim, we used some a priori knowledge about the goal direction, defined between x 0 and x g , provided by the isotropic FMM using Vesselness(x) cost function (see Fig. 3(a)). Indeed, this FMM step provides the distance map U c between x 0 and x g (or reciprocally), and then the goal direction field ∇U c (x). This helps to define the appropriate guidance field G(x) ∝ ∇U c (x). Some variations could be added on the goal direction field to enlarge the guidance field set. The admissible control set A can then be defined as:

A = {A × G}, (10)
where A ⊂]0; a max ] is a finite set of action magnitudes, and G = {G(x) ∼ [0; π 2 ]} is the guidance field set. Let us recall that A could be filtered to consider some technological constraints (such as the duty-cycle, or the magnetic slew rate).

This planning framework has been validated on the case study dataset considering different microrobot sizes navigating within arteries. First, Fig. 8 demonstrates the influence of the motion force magnitude on finding the minimum-feasibleenergy path. We noticed that as the motion force magnitude increased, the required energy E ϖ (P) increased exponentially. This is significant for important microrobot size (r = 750 µm on Fig. 8(a)). Indeed, Fig. 8(b) shows that with much power the microrobot is more able to leave the vessel's centerline P c (provided by the isotropic FMM) to reach the goal faster. If limited magnetic gradient magnitudes are available to propel the microrobot, such as in clinical MRI, Fig. 9 shows that increasing its size should overcome the power limitation. Therefore, the proposed min-energy path planning has been evaluated for various microrobot sizes and for EMA and MRN systems, as reported on Fig. 10. We can notice that for a microrobot radius bigger than 1.5 mm the curve length increases slightly, and starts to decrease significantly below 1 mm (cf. Fig. 10(a)). Furthermore, Fig. 10(b) depicts that the energy functional decreases significantly with the microrobot radius. In fact, the action force f a (a) decreases following a cubic power (L 3 ), minimizing the power budget. Therefore, it is more difficult for neodymium microrobots with a size below 750 µm to overcome the blood flow constraints that take place in the artery. Finally, as an illustration of the blood flow influence, Fig. 11 shows the planning results of our minimumfeasible-path finder for three targeted points (x g1 , x g2 and x g3 ), and Table I lists the curve lengths of their geodesic (P 1 , P 2 and P 3 ) and the energy functional E ϖ (P). The min-energy geodesic is curve P 3 , while the shortest length is P 2 . This means that although P 3 is longer than P 2 , P 3 consumes less energy.

IV. DISCUSSION

The most computationally intensive part of the microrobotic surgical application is the vesselness filter (that provides Vesselness(x) cost function) based on a multiscale analysis of the Hessian matrix [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF]. This multiscale interpretation of eigenvalues of the Hessian matrix could be improved thanks to parallelization through GPU implementation. Furthermore, as mentioned, instead of using the overall dataset grid (in the case study N = 147 456), the vector fields simulations are performed on only 18 941 nodes corresponding to the vasculature. Thus, on an Intel Core2 3Ghz processor, one simulation time takes less than 150 ms. For the periodic pulsatile blood flow, this allows to compute one cardiac cycle with 50 time-steps in less than 8 s. Therefore, the combination of the microrobot and the caradiovascular model together with the FMM may become accurate enough to support interactive microrobotic surgical planning applications. Finally, in the planning step, the main computational aspect remains on action set A generation. This set A is related to the dimension of the action magnitude set A and guidance field set G, specified by the users. For the case study with dim A = {10 × 10}, the overall planning procedure needs a calculation time of about 36 s. This computational efficiency ensures a high level of interactivity for preoperative simulation, planning and navigation.

In the simulation process we have chosen to focus on computational efficiency, rather than on accurate blood flow simulations. Nonetheless, the vasculature and blood flow modeling are still strongly investigated by the communities, and remain an open issue. In future extends, we will use some analytical modeling of vessel bifurcation as in [START_REF] Arcèse | Endovascular Magnetically-Guided Robots: Navigation Modeling and Optimization[END_REF], to improve the blood flow simulation accuracy. Although the visualization of vector fields is still a challenging issue (especially in 3D), more intuitive flow rendering (such as streamline, streamribbon, etc.) could also be helpful for the microrobotic surgery applications.

V. CONCLUSION In this study, a microrobotic surgical simulation and planning has been presented for the treatment of peripheral arterial diseases. However, the proposed framework could also be applied to other regions of the body, different therapy and various MIS microrobots. Especially, the results shows that interactive applications can be helpful to design the suitable microrobot design (e.g. shape, size, materials, etc.) wrt. anatomic consistency to improve microrobotic surgical applications. Hence, preoperative surgical planning with blood flow simulation provides useful information for understanding cardiovascular pathologies, predicting their onset and choosing an optimal therapy strategy for innovative microrobotic MIS procedure.

Fig. 1 .

 1 Fig. 1. Modeling scheme example in case of bead pulled microrobot.

Fig. 2 .

 2 Fig. 2. Microforces evolution: (a) apparent weight fg (b) magnetic force fa (with ∇B = 200 mT/m), (c) drag force f d (with ∥v f ∥ = 400 mm/s), (d) electrostatic force fe, (e) van der Waals force fv, as function of the distance to the vessel wall δ in artery (R = 5 mm) and wrt. the microrobot radius r. (f) Microforces acting on a microrobot with r = 1 mm wrt. the distance to the vessel wall δ in artery (R = 5 mm).

Fig. 3 .

 3 Fig. 3. Preoperative microrobot surgical simulation and planning application: (a) application workflow; (b) vessel centerline extraction example.

Fig. 4 .

 4 Fig. 4. Real MIP image from MRA (resolution of N = 384 × 384 nodes) obtained in a patient with an isolated stenosis in the iliac artery (the white (A)arrow), used to illustrate microrobotic MIS treatment simulation and planning.

Fig. 5 .

 5 Fig. 5. Simulation results of the blood flow velocity field computation.

Fig. 6 .

 6 Fig. 6. Simulation results of the force field: f = f d (v) + fe(x) + fg, applied on a microrobot of r = 500 µm radius.

Fig. 7 .

 7 Fig. 7. Pulsatile blood flow simulation: (a) the periodic time-varying velocity in the descending aorta; (b) zoomed view of the blood flow velocity field; (cd) zoomed view of the force field f = f d + fe + fg acting on a microrobot of: (c) r = 250 µm, and (d) r = 1500 µm for different time instants.

Fig. 8 .

 8 Fig. 8. Planning results on varying control action influence for a neodymium microrobot of size 250 µm and 750 µm: (a) the energy functional Eϖ(P) wrt. the action magnitude (i.e. the magnetic gradient) and (b) zoomed view of the corresponding geodesics: Pc the vessel centerline and Pa the min-energy path for the varying action magnitude.

Fig. 9 .

 9 Fig. 9. Planning results on varying control action influence for a neodymium microrobot of size 1.5 mm and 2.5 mm in clinical MRI context.

a

  Fig. 10. Planning results with varying microrobot radii, for an EMA system (amax = 350 mT/m) and a MRN (amax = 80 mT/m): (a) the geodesic curvature length (in mm); and (b) the energy functional.

Fig. 11 .

 11 Fig. 11. Planning results using the proposed anisotropic cost function for a neodymium microrobot with r = 1.5 mm: three targeted point are specified (x g1 , x g2 and x g3 ). The colored map depicted the power functional level.

TABLE I PLANNING

 I

	Geodesic #	Length (mm)	Energy Eϖ(P)
	P 1	360.086	0.6765
	P 2	359.612	0.5118
	P 3	361.735	0.4879

RESULTS OF MINIMUM-FEASIBLE-PATH FOR THREE TARGETED POINT (x g1 , x g2 AND x g3 ) DEPICTED ON FIG.11

ACKNOWLEDGMENT

The PhD student K. Belharet is supported by the Nano-IRM project, founded by Region Centre and City of Bourges.