Samir Adly
email: samir.adly@unilim.fr

Huynh Van Ngai

Nguyen Van Vu
email: van-vu.nguyen@etu.unilim.fr

Antonino Maugeri

Newton-Type Method for Solving Generalized Equations on Riemannian Manifolds

Keywords: Riemannian manifold, generalized equation, Newton's method, metric regularity, variational inclusion

Dedicated to

Introduction

We consider the problem finding a solution of inclusion 0 ∈ f (x) + F (x).

(1)

Here, the variable x varies in a finite dimensional Riemannian manifold M, f : M -→ R n is at least continuous, while F : M ⇒ R n is a set-valued mapping. We suppose that the graph Gr(F ) of F is closed with respect to the product topology on M × R n .

The general inclusion (1) covers many situations which have been studied widely in the literature. When M is a Euclidean space, (1) is nothing but the so-called generalized equation. In [START_REF] Adly | Newton's method for solving inclusions using set-valued approximations[END_REF], the authors studied the (super-)linear convergence of a Newton-type iterative process for solving generalized equations. The Kantorovich approach and Smale's classical (α, γ)-theory was extended to generalized equations in [START_REF] Adly | Newton's method for solving generalized equation: Kantorovich's and Smale's approaches[END_REF]. If F (x) ≡ 0, problem (1) reduces to solve nonlinear equation f (x) = 0 on M. In the case F (x) ≡ K for a fixed cone K ⊂ R n , (1) coincides with the problem studied in [START_REF] Wang | Extended Newton's method for mappings on Riemannian manifolds with values in a cone[END_REF].

The current work considers a scheme of Newton-type method to approximate a solution of [START_REF] Absil | Sepulchre: Optimization Algorithms on Matrix Manifolds[END_REF]. This is based on the well-known Josephy-Newton method applied to the generalized equation which was introduced in [START_REF] Dontchev | Implicit Functions and Solution Mappings: A View from Variational Analysis[END_REF]. The strategy is that to start at a guess point x 0 nearby a solution, and generates a sequence (x k , v k ) in the tangent bundle T M by the scheme

0 ∈ f (x k ) + Df (x k )(v k ) + (F • R x k ) (v k ), x k+1 = R x (v k ). (2)
In [START_REF] Adler | Newton's method on Riemannian manifolds and a geometric model for the human spine[END_REF], Df (x) : T x M -→ R n and R x : T x M -→ M are respectively the covariant derivative of f and the retraction at x (see Section 2). Observe that when F ≡ 0, (2) subsumes as a particular case of Newton-type methods for solving nonlinear equation on manifold M studied e.g. in [START_REF] Adler | Newton's method on Riemannian manifolds and a geometric model for the human spine[END_REF][START_REF] Dedieu | Newton's method on Riemannian manifolds: covariant alpha theory[END_REF]. Furthermore, it might be viewed as an extension of algorithms for finding singularities of a vector field considered by the works [START_REF] Alvarez | A unifying local convergence result for Newton's method in Riemannian manifolds[END_REF][START_REF] Ferreira | Local convergence of Newton's method under a majorant condition in Riemannian manifolds[END_REF][START_REF] Ferreira | Kantorovich's theorem on Newton's method in Riemannian manifolds[END_REF] if f is replaced by a smooth vector field X, and F is the zero field F (x) = 0 x . In the same spirit, more other discussions about Newton-type method applying on smooth manifold can be found in [START_REF] Absil | Sepulchre: Optimization Algorithms on Matrix Manifolds[END_REF][START_REF] Smith | Optimization techniques on Riemannian manifolds, in: Hamiltonian and Gradient Flows[END_REF][START_REF] Ring | Optimization methods on Riemannian manifolds and their application to shape space[END_REF][START_REF] Udrişte | Convex Functions and Optimization Methods on Riemannian Manifolds[END_REF].

The paper is organized as follows. In section 2, we recall basic elements, notations and backgrounds on Riemannian geometry that will be useful in the sequel. Section 3 is devoted to the stability of the metric regularity property of the sum of two operators on Riemannian manifolds. In section 4, we prove the local and global convergence of the retraction Newton-type algorithm for solving [START_REF] Absil | Sepulchre: Optimization Algorithms on Matrix Manifolds[END_REF].

Notions and backgrounds

Throughout this paper, we prefer to adopt the standard notations and different concepts used in [START_REF] Carmo | Riemannian Geometry, Mathematics: Theory and Applications[END_REF][START_REF] Lee | Riemannian manifolds: An Introduction to Curvature[END_REF] dealing with the basic background from Riemannian geometry. For more details about these events, the readers are referred to [START_REF] Helgason | Differential Geometry, Lie Groups, and Symmetric Spaces[END_REF][START_REF] Absil | Sepulchre: Optimization Algorithms on Matrix Manifolds[END_REF][START_REF] Lang | Fundamentals of Differential Geometry[END_REF][START_REF] Sakai | Riemannian Geometry[END_REF] and references therein.

Riemannian manifold, metric structure. All objects under consideration in this paper concerns finite dimensional smooth manifolds. The terminology smooth is meant to be differentiable at least beyond the order that appears. Elementary notions of differential geometry are assumed to be familiar.

A Riemannian manifold M of dimension m is a differentiable manifold of the same dimension which is endowed with a Riemannian metric g. If x is a point in M, the norm induced by g in the tangent space T x M of M at x is denoted by ∥•∥ x . By B x (v, r) (resp. B x (v, r)) we mean the open (resp. closed) ball in T x M with center at v ∈ T x M and radius r > 0. To mention about the open (closed) unit ball in T x M, we write B x (resp. B x ).

Let χ : [a, b] -→ M be a piecewise smooth curve, then its length is given by the quantity

ℓ(χ) = ∫ b a ∥χ ′ (t)∥ χ(t) dt.
Recall that here χ ′ (t) stands for the tangent vector (or velocity) of χ at the instant t. (This is sometimes denoted in other ways, such as, χ(t) or dχ dt (t).) For two points x, y in the manifold M, the Riemannian distance between x and y can be defined as follow

d M (x, y) = inf { ℓ(χ) χ : [a, b] ⊂ R -→ M piecewise smooth; χ(a) = x, χ(b) = y } .
We will omit the subscript M when the manifold is fixed.

With the Riemannian distance d M , M becomes a metric space. Unless it has some other specification, the space (M, d M ) is always supposed to be complete.

Completeness can be described via the famous Hofp-Rinow theorem (see [START_REF] Carmo | Riemannian Geometry, Mathematics: Theory and Applications[END_REF][START_REF] Lee | Riemannian manifolds: An Introduction to Curvature[END_REF]). As usual, we denote the open and closed ball on M with center x and radius r > 0 by B M (x, r) and B M (x, r), respectively.

Connection, covariant derivative. Let ∇ be the Levi-Civita connection on the manifold M. Consider a smooth function f : M -→ R and a vector field Y on M. The covariant derivative of f with respect to Y is the function

∇ Y f := Y f
, where Y f indicates the action of Y on f (see [START_REF] Helgason | Differential Geometry, Lie Groups, and Symmetric Spaces[END_REF][START_REF] Lee | Riemannian manifolds: An Introduction to Curvature[END_REF]). For x ∈ M, the covariant derivative Df (x) :

T x M -→ R of f at x is a function which assigns to each v ∈ T x M the value Df (x)(v) := (∇ Y f )(x) ∈ R, where Y is a vector field such that Y x = v. Explicitly, we have Df (x)(v) = χ ′ (0)(f ) = (f • χ) ′ (0) (3)
for any smooth curve χ satisfying χ(0) = x and χ ′ (0

) = v. Generally, if f = (f 1 , . . . , f n ) : M -→ R n is a smooth map, then one defines Df (x) =    Df 1 (x) . . . Df n (x)   
as the covariant derivative of f at x. The covariant derivative Df (x) is a linear map from T x M into R n (cf. [START_REF] Dedieu | Newton's method on Riemannian manifolds: covariant alpha theory[END_REF]). In the sequel, we use the following norm for the covariant derivative

∥Df (x)∥ := sup u∈TxM ∥u∥ x 1 ∥Df (x)(u)∥ R n . (4)
Vector transportation. We will need later the concept of vector transports. They provide a way to link between different tangent spaces of a manifold. Among them, parallel transportations are the most typical and important. Assume χ : [a, b] -→ M is a smooth curve. By a vector field V along χ, we mean a smooth map such that V (t) ∈ T χ(t) M. V is said to be parallel along χ if its covariant derivative ∇ χ ′ V vanishes. Then, the parallel transport

P χ,χ(a),χ(b) : T χ(a) M -→ T χ(b) M is given by P χ,χ(a),χ(b) (v) = V (b),
where V is the unique vector field along χ satisfying

∇ χ ′ V = 0 and V (a) = v.
It is well-known from Riemannian geometry that P χ,χ(a),χ(b) is a linear isometry.

In particular, one always has

P χ,χ(a),χ(b) T χ(a) M,T χ(b) M = 1 and P -1 χ,χ(a),χ(b) = P χ,χ(b),χ(a) .
Retraction. Retraction (cf. [START_REF] Alvarez | A unifying local convergence result for Newton's method in Riemannian manifolds[END_REF][START_REF] Ring | Optimization methods on Riemannian manifolds and their application to shape space[END_REF]) is a crucial object for our approach. In this paper, a retraction is meant to be a smooth map R : T M -→ M from the tangent bundle into M so that its restriction R x (retraction at x) of R onto each tangent space

T x M satisfies • R x (0 x ) = x, • (dR x ) 0x = id TxM , with the identification T 0x (T x M) ≃ T x M.
On the preceding descriptions, (dR x ) 0x stands for the differential of R x at the origin 0 x of T x M. If R is a retraction, and (x, v) is in T M, then the map χ : t -→ R x (tv) defines a smooth curve on M satisfying χ(0) = x and χ ′ (0) = v. There is a natural retraction reduced by geodesics on M, known as exponential map of the tangent bundle [START_REF] Carmo | Riemannian Geometry, Mathematics: Theory and Applications[END_REF]. Recall that a curve

χ : I -→ M is called a geodesic, if its acceleration ∇ χ ′ χ ′ is vanishing. The exponential map exp : T M -→ M can be defined by exp(x, v) = χ(1, x, v), where t -→ χ(t, x, v) is the unique geodesic such that χ(0, x, v) = x and χ ′ (0, x, v) = v.
The exponential map has many important properties. At each point of a Riemannian manifold, there is a normal neighborhood for which exp x is injective on some ball around the origin of tangent space T x M (see [START_REF] Helgason | Differential Geometry, Lie Groups, and Symmetric Spaces[END_REF][START_REF] Carmo | Riemannian Geometry, Mathematics: Theory and Applications[END_REF][START_REF] Sakai | Riemannian Geometry[END_REF]). In accordance with this property, we will establish here a similar result applied for any retraction. The proof is based on the same strategy as the case of exponential maps presented in the literature.

Proposition 2.1. Let M be a complete m-dimensional Riemannian manifold, and let R : T M -→ M be a retraction. Then for every x ∈ M there exists a neighborhood W x = B M (x, r x ) (r x > 0 depends on x) and a real number

C R (x) > 0 such that if y ∈ W x then R y : T y M -→ M is injective in the ball C R (x)B y with R y (C R (x)B y ) ⊃ W x .
We will call such a pair (r x , C R (x)) as a R-normal pair at x. Proof. Let (U, x) be a local coordinate of M at x. Then

( U × U, (x, x)
) is a local coordinate on M × M. Consider the map H : T U -→ M × M given by

H(z, v) = ( z, R z (v) ) . Since R x (0 x ) =
x and (dR x ) 0x = id, the matrix of dH (x,0x) in the local coordinate above can be written as follow

( I 0 * I
) .

Hence, we can apply an analogous argument as the proof of existence for normal neighborhood in Riemannian manifold and obtain the required conclusion. For shortness, one keeps in mind [15, Lemma 5.12] and/or [START_REF] Carmo | Riemannian Geometry, Mathematics: Theory and Applications[END_REF]Theorem 3.7].

Besides Proposition 2.1, we shall also need some other facts. The next statement will be useful throughout the rest of this paper. M and positive numbers ε, ρ 1 , ρ 2 . We say that R satisfies the uniform rate condition on Ω with respect to ratios ρ 1 , ρ 2 and radius ε, written as R ∈

URC(ρ 1 , ρ 2 , ε, Ω), if ρ 1 ∥u -v∥ x d ( R x (u), R x (v) ) ρ 2 ∥u -v∥ x (5)
for all x ∈ Ω and u, v ∈ εB x .

Remark 2.3. It is obvious to see that [START_REF] Alvarez | A unifying local convergence result for Newton's method in Riemannian manifolds[END_REF] is trivial in the case M = R n and R is the usual translation R x (u) = x + u. Another less trivial case will be showed in Example 2.4 below. When R = exp is the exponential, the left-hand side of ( 5) is satisfied on any Hadamard manifold (for instance, cf. [START_REF] Sakai | Riemannian Geometry[END_REF], Chapter V, Proposition 4.5). If R = exp and Ω = B M (x, r), where r is small enough, the right-hand side of ( 5) is also valid in terms of local property [START_REF] Helgason | Differential Geometry, Lie Groups, and Symmetric Spaces[END_REF][START_REF] Dedieu | Newton's method on Riemannian manifolds: covariant alpha theory[END_REF]. For general retractions, a condition ensures (5) appeared in the work [START_REF] Ring | Optimization methods on Riemannian manifolds and their application to shape space[END_REF], where retractions R x are required to satisfy the equicontinuous derivative on Ω.

Example 2.4 (retraction on the unit sphere). Consider the m-unit sphere

S m := { x ∈ R m+1 : ∥x∥ 2 = x T x = 1 }
, where the vector is written in column form. S m is endowed with the Riemannian metric

g x (u, v) = ⟨u, v⟩ TxS m := u T v, ∀u, v ∈ T x S m = { z ∈ R m+1 : z T x = 0 } . (6)
The Riemannian distance associated to the metric above is given by

d(x, y) = arccos ( x T y ) , x, y ∈ S m . (7)
According to [START_REF] Absil | Sepulchre: Optimization Algorithms on Matrix Manifolds[END_REF], let us consider the retraction

R x (u) = x + u ∥x + u∥ , x ∈ S m , u ∈ T x S m . (8)
Observe that we have

d ( x, R x (u) ) = arccos ( x T (x + u) ∥x + u∥ ) = arccos 1 √ 1 + ∥u∥ 2 = arctan ∥u∥ .
So, R x is injective on the whole space T x S m and R x (εB x ) = B S m (x, arctan ε) for every x ∈ S m and ε > 0.

Let r > 0 be small enough and Ω ⊂ S m is a closed set. Suppose that x ∈ Ω and u, v ∈ rB x ⊂ T x S m are fixed. For brevity, we put

d := d ( R x (u), R x (v)
) . From ( 7) and ( 8) we get

cos d = R x (u) T R x (v) = 1 ∥x + u∥ ∥x + v∥ (x + u) T (x + v) = 1 + u T v ( 1 + ∥u∥ 2 ) 1/2 ( 1 + ∥v∥ 2 ) 1/2 .
Hence

sin 2 d = 1 -cos 2 d = ∥u -v∥ 2 + ∥u∥ 2 ∥v∥ 2 - ( u T v ) 2 ( 1 + ∥u∥ 2 ) ( 1 + ∥v∥ 2 ) . It is clear that ∥u∥ 2 ∥v∥ 2 - ( u T v ) 2 0, so sin 2 d ∥u -v∥ 2 ( 1 + ∥u∥ 2 ) ( 1 + ∥v∥ 2 ) ( 1 1 + r 2 ) 2 ∥u -v∥ 2 . (9)
On the other hand, we have 2u

T v = ∥u∥ 2 + ∥v∥ 2 -∥u -v∥ 2 , which implies ∥u∥ 2 ∥v∥ 2 - ( u T v ) 2 = ∥u∥ 2 ∥v∥ 2 - 1 4
( ∥u∥ 2 + ∥v∥ 2 -∥u -v∥ 2 ) 2 = 1 2
( ∥u∥ 2 + ∥v∥ 2 ) ∥u -v∥ 2 - 1 4 [ ∥u -v∥ 4 + ( ∥u∥ 2 -∥v∥ 2 ) 2 ]
1 2

( ∥u∥ 2 + ∥v∥ 2 ) ∥u -v∥ 2 .
Thus,

sin 2 d [ 1 + 1 2 ( ∥u∥ 2 + ∥v∥ 2 )] ∥u -v∥ 2 ( 1 + ∥u∥ 2 ) ( 1 + ∥v∥ 2 ) ∥u -v∥ 2 . ( 10
)
Note that

d d ( x, R x (u) ) + d ( x, R x (v) )
= arctan ∥u∥ + arctan ∥v∥ π 2 whenever u, v are in the unit ball of T x S m . Taking into account the following fact

0 < t π 2 =⇒ 2 π sin t t 1, (9)
and [START_REF] Ferreira | Kantorovich's theorem on Newton's method in Riemannian manifolds[END_REF] give us

1 1 + r 2 d ( R x (u), R x (v) ) ∥u -v∥ π 2 , 0 r 1, u ̸ = v ∈ rB x .
Consequently,

R ∈ URC (ρ 1 (r), ρ 2 (r), r, Ω) , (11)
where r ∈ (0, 1], ρ 1 (r) = 1 1+r 2 and ρ 2 (r) = π 2 .

Stability of the metric regularity property on Riemannian manifolds

To study the convergence of schemes like (2), the concept of metric regularity for set-valued mapping [START_REF] Ioffe | Metric regularity and subdifferential calculus[END_REF] plays an important role in our analysis. Recall that a set-valued (or also multivalued) mapping Φ : X ⇒ Y between two metric spaces X and Y is a correspondence assigning to each x ∈ X a subset Φ(x) ⊂ Y . Such a mapping is totally determined by its graph Gr(Φ) :=

{ (x, y) ∈ X × Y : y ∈ Φ(x)
} . The mapping Φ is said to be metrically regular on a subset V ⊂ X × Y if there exists a constant κ > 0 (a modulus of regularity) such that

dist X ( x, Φ -1 (y) ) κdist Y ( y, Φ(x) ) , for all (x, y) ∈ V. (12)
In [START_REF] Helgason | Differential Geometry, Lie Groups, and Symmetric Spaces[END_REF], dist X (•, •) and dist Y (•, •) are respectively the distances in X and Y , while Φ -1 stands for the inverse of Φ, given by

x ∈ Φ -1 (y) ⇐⇒ y ∈ Φ(x).
When dealing with a property like ( 12), we will sometimes write κ ∈ REG ( Φ, V ) to simplify notations. Moreover, we also write G ∈ L ip λ (S) to indicate that the mapping (maybe set-valued) G is Lipschitz continuous on the set S with a modulus λ. In words,

G ∈ L ip λ (S) means d H Y ( G(x), G(y) ) λd X (x, y), x, y ∈ S, (13)
where

d H Y ( C, D ) is the Hausdorff distance [16] between two subsets C and D in Y d H Y ( C, D ) = sup y∈Y d Y ( y, C ) -d Y ( y, D ) . ( 14
)
We now present the main results of this section, which ensures the stability of metric regularity property. These are essential preliminaries for the convergence analysis at the end of this paper.

Proposition 3.1 (stability of local metric regularity). Given a complete

Riemannian manifold of dimension m, and a retraction R : T M -→ M. Fix a point x ∈ M and let (r x , C R (x)) be a R-normal pair at x. Suppose that there exist some positive numbers 0

< ρ 1 ρ 2 satisfying R ∈ URC(ρ 1 , ρ 2 , δ, W x ) with W x = B M (x, r x ) and δ = C R (x)
. Pick some tangent vector ū ∈ T x M and some positive constants r, s, r ′ , s ′ , λ, λ ′ , σ and κ obeying the following relations

                   θ = ρ -1 1 ρ 2 κ (λ ′ + λµ x ) < 1, ρ -1 1 ρ 2 ( 1 + 1+θ 1-θ ) r ′ + κ 1-θ s ′ r, ( ρ -1 1 ρ 2 ) 2 (λ ′ + λµ x ) ( 1 + 1+θ 1-θ ) r ′ + ( 1 + θ 1-θ ) s ′ s, ρ -1 1 ρ 2 ∥ū∥ x + ρ -1 1 σ + ρ -1 1 ρ 2 ( 1 + 1+θ 1-θ ) r ′ + κ 1-θ s ′ δ, ρ 2 ∥ū∥ x + ρ 2 [ ρ -1 1 ρ 2 ( 1 + 1+θ 1-θ ) r ′ + κ 1-θ s ′ ] r x , (15)
for

µ x = ∥Df (x)∥. Assume that the mapping Φ x = Df (x)(•) + (F • R x ) (•) is metrically regular on a neighborhood V x := B x (ū, r) × B R n (p, s) of (ū, p) ∈ Gr(Φ x ) together with a modulus κ. Pick some point y ∈ B M (x, σ) ∩ W x such that, for each geodesic χ : [0, 1] -→ M having χ(0) = x, χ(1) = y and χ([0, 1]) ⊂ W x : (i) Σ χ,y,x = R -1 x • R y -P χ,y,x ∈ L ip λ (δB y ); (ii) the linear map G χ,x,y := Df (y) • P χ,x,y -Df (x) satisfies ∥G χ,x,y ∥ λ ′ . Set v = ( R -1 y • R x ) (ū) ∈ T y M, q = p -Df (x)(ū) + Df (y)(v) ∈ R n , and τ = ρ -1 1 ρ 2 κ 1-θ . Then one has τ ∈ REG ( Φ y , V y ) , where Φ y = Df (y)(•) + (F • R y ) (•) and V y := B y (v, r ′ ) × B R n (q, s ′ ).
To prove Proposition 3.1, we need the following lemma.

Lemma 3.2. Under the assumptions of Proposition 3.1, one has

Ψ y,x := Df (y) -Df (x) • ( R -1 x • R y ) ∈ L ip λ ′ +λµx (δB y ).
Proof. Without loss of generality, we can focus on the case where W x is a normal neighborhood at x. Let χ : [0, 1] -→ W x be a geodesic having χ(0) = x and χ(1) = y. It is not difficult to check that

Ψ y,x (v) = ( G χ,x,y • P -1 χ,x,y ) (v) + (Df (x) • Σ χ,y,x ) (v),
for v ∈ δB y . As a result, we find

∥Ψ y,x (v ′ ) -Ψ y,x (v)∥ R n = ( G χ,y,x • P -1 χ,x,y ) (v ′ -v) + Df (x) (Σ χ,y,x (v ′ ) -Σ χ,y,x (v)) R n ∥G χ,y,x ∥ P -1 χ,x,y ∥v ′ -v∥ y + ∥Df (x)∥ ∥Σ χ,y,x (v ′ ) -Σ χ,y,x (v)∥ x λ ′ ∥v ′ -v∥ y + µ x λ ∥v ′ -v∥ y .
Hence, the conclusion of Lemma 3.2 follows.

Proof of Proposition 3.1. Pick (v, w) ∈ V y . We have to establish the estimation

dist ( v, Φ -1 y (w) ) τ dist ( w, Φ y (v) ) . ( 16
)
We omitted the subscripts of spaces in [START_REF] Papadopoulos | Metric Spaces, Convexity and Nonpositive Curvature[END_REF], since these ones are determined by the objects upon which the corresponding distances act. Here and in what follows, we will use the common notation dist(•, •) for distance on any arbitrary metric space. The set Φ y (v) is obviously closed. So, ( 16) is easy when dist ( w, Φ y (v)

) = 0. Consider the case η = dist ( w, Φ y (v) ) > 0. First, we set v 0 = v and y 0 = R y (v 0 ). Note that v = R -1 y (R x (ū)) ∈ C R (x)B y ,

the triangle inequality and Assumption 2.2 give us

∥v∥ y ρ -1 1 d ( y, R y (v) ) = ρ -1 1 d ( y, R x (ū) ) ρ -1 1 ( d(y, x) + d ( x, R x (ū) )) ρ -1 1 (d(x, y) + ρ 2 ∥ū∥ x ) ρ -1 1 σ + ρ -1 1 ρ 2 ∥ū∥ x . This implies ∥v 0 ∥ y ∥v 0 -v∥ y + ∥v∥ y < r ′ + ρ -1 1 σ + ρ -1 1 ρ 2 ∥ū∥ x < δ. Since R ∈ URC (ρ 1 , ρ 2 , δ, W x ), we have d ( y 0 , x ) d ( y 0 , R y (v) ) + d ( R y (v), x ) = d ( R y (v 0 ), R y (v) ) + d ( R x (ū), x ) ρ 2 ∥v 0 -v∥ y + ρ 2 ∥ū∥ x < ρ 2 r ′ + ρ 2 ∥ū∥ x < ρ 2 ∥ū∥ x + ρ 2 [ ρ -1 1 ρ 2 ( 1 + 1 + θ 1 -θ ) r ′ + κ 1 -θ s ′ ] r x .
Thus,

y 0 = R y (v 0 ) ∈ W x ⊂ R x (δB x ).
Let us choose u 0 ∈ δB x with R x (u 0 ) = y 0 and put

w 0 = w + Df (x)(u 0 ) -Df (y)(v 0 ) ∈ R n .
We claim that (u 0 , w 0 ) ∈ V x . Indeed, due to Assumption 2.2 and the choice of

u 0 we get ∥u 0 -ū∥ x ρ -1 1 d ( R x (u 0 ), R x (ū) ) = ρ -1 1 d ( R y (v 0 ), R y (v) ) ρ -1 1 ρ 2 ∥v 0 -v∥ y , (17)
which yields ∥u 0 -ū∥

x < ρ -1 1 ρ 2 r ′ < r. Taking into account R x (u 0 ) = R y (v 0 ) and R x (ū) = R y (v), we deduce ∥w 0 -p∥ R n = ∥w -Ψ y,x (v 0 ) -[q -Ψ y,x (v)]∥ R n ∥Ψ y,x (v 0 ) -Ψ y,x (v)∥ R n + ∥w -q∥ R n (λ ′ + λµ x ) ∥v 0 -v∥ + ∥w -q∥ R n . (18)
In view of [START_REF] Sakai | Riemannian Geometry[END_REF], we obtain

∥w 0 -p∥ R n < (λ ′ + λµ x )r ′ + s ′ < ( ρ -1 1 ρ 2 ) 2 (λ ′ + λµ x ) ( 1 + 1 + θ 1 -θ ) r ′ + ( 1 + θ 1 -θ ) s ′ s.
Hence, the inclusion

(u 0 , v 0 ) ∈ V x is clear. By invoking the fact κ ∈ REG(Φ x , V x ), we find dist ( u 0 , Φ -1 x (w 0 ) ) κdist ( w 0 , Φ x (u 0 ) ) = κdist ( w 0 , Df (x)(u 0 ) + (F • R x ) (u 0 ) ) = κdist ( w, Df (y)(v 0 ) + (F • R y ) (v 0 ) ) = κdist ( w, Df (y)(v) + (F • R y ) (v) ) = κη.
So we can select a vector u 1 in Φ -1 x (w 0 ) which satisfies

∥u 0 -u 1 ∥ x = dist ( u 0 , Φ -1 x (w 0 ) ) κη.
To continue, let us set y 1 = R x (u 1 ). Inasmuch as κ ∈ REG(Φ x , V x ), one has

dist ( ū, Φ -1 x (w 0 ) ) κdist ( w 0 , Φ x (ū) ) κ ∥w 0 -p∥ R n .
Consequently,

∥u 1 -u 0 ∥ x = dist ( u 0 , Φ -1 x (w 0 ) ) ∥u 0 -ū∥ x + dist ( ū, Φ -1 x (w 0 ) ) ∥u 0 -ū∥ x + κ ∥w 0 -p∥ R n . (19)
By induction hypothesis, suppose that the tangent vectors u 0 , . . . , u k ∈ B x (ū, r), v 0 , . . . , v k-1 ∈ δB y ⊂ T y M are given, and they obey the following conditions

• R y (v j ) = R x (u j ), j = 0, 1, . . . , k -1; • u j+1 ∈ Φ -1 x (w j ), with w j = w+Df (x)(u j )-Df (y)(v j ), j = 0, 1, . . . , k-1; • ∥u j -u j+1 ∥ x θ j ∥u 1 -u 0 ∥ x , for θ = ρ -1 1 ρ 2 κ(λ ′ + λµ x ) < 1 and j = 0, 1, . . . , k -1.
Thanks to the triangle inequality, we can write

∥u j -ū∥ x j-1 ∑ i=0 ∥u i -u i+1 ∥ x + ∥u 0 -ū∥ x j-1 ∑ i=0 θ i ∥u 1 -u 0 ∥ x + ∥u 0 -ū∥ x = 1 -θ j 1 -θ ∥u 1 -u 0 ∥ x + ∥u 0 -ū∥ x 1 1 -θ ∥u 1 -u 0 ∥ x + ρ -1 1 ρ 2 ∥v 0 -v∥ y .
Since ρ -1 1 ρ 2 1, from ( 17), ( 18) and ( 19) we obtain

∥u j -ū∥ x { ρ -1 1 ρ 2 + 1 1 -θ [ ρ -1 1 ρ 2 + κ(λ ′ + λµ x ) ] } ∥v 0 -v∥ x + κ 1 -θ ∥w -q∥ R n ρ -1 1 ρ 2 ( 1 + 1 + θ 1 -θ ) ∥v 0 -v∥ x + κ 1 -θ ∥w -q∥ R n .
The latter implies

∥u j ∥ x ∥u j -ū∥ x + ∥ū∥ x ρ -1 1 ρ 2 ( 1 + 1 + θ 1 -θ ) ∥v 0 -v∥ x + κ 1 -θ ∥w -q∥ R n + ∥ū∥ x < ρ -1 1 ρ 2 ( 1 + 1 + θ 1 -θ ) r ′ + κ 1 -θ s ′ + ∥ū∥ x δ.
Next, involving again Assumption 2.2, we derive

d ( R x (u k ), x ) d ( R x (u k ), R x (ū) ) + d ( R x (ū), x ) ρ 2 ∥u k -ū∥ x + ρ 2 ∥ū∥ x ρ 2 [ ρ -1 1 ρ 2 ( 1 + 1 + θ 1 -θ ) ∥v 0 -v∥ x + κ 1 -θ ∥w -q∥ R n ] + ρ 2 ∥ū∥ x < ρ 2 [ ρ -1 1 ρ 2 ( 1 + 1 + θ 1 -θ ) r ′ + κ 1 -θ s ′ ] + ρ 2 ∥ū∥ x r x .
This means

y k = R x (u k ) ∈ W x . Selecting v k ∈ δB y with R y (v k ) = R x (u k )
and

w k = w + Df (x)(u k ) -Df (y)(v k ) = w -Ψ y,x (v k ),
we get

∥w k -p∥ R n = ∥w -Ψ y,x (v k ) -[q -Ψ y,x (v)]∥ R n ∥Ψ y,x (v k ) -Ψ y,x (v)∥ R n + ∥w -q∥ R n (λ ′ + λµ x ) ∥v k -v∥ y + ∥w -q∥ R n < (λ ′ + λµ x ) ∥v k -v∥ y + s ′ .
Moreover, it holds that

∥v k -v∥ y ρ -1 1 d ( R y (v k ), R y (v) ) = ρ -1 1 d ( R x (u k ), R x (ū) ) ρ -1 1 ρ 2 ∥u k -ū∥ x ρ -1 1 ρ 2 { ρ -1 1 ρ 2 ( 1 + 1 + θ 1 -θ ) ∥v 0 -v∥ x + κ 1 -θ ∥w -q∥ R n } < ρ -1 1 ρ 2 { ρ -1 1 ρ 2 ( 1 + 1 + θ 1 -θ ) r ′ + κ 1 -θ s ′ } .
Hence,

∥w k -p∥ R n < ρ -1 1 ρ 2 (λ ′ + λµ x ) { ρ -1 1 ρ 2 ( 1 + 1 + θ 1 -θ ) r ′ + κ 1 -θ s ′ } + s ′ = ( ρ -1 1 ρ 2 ) 2 (λ ′ + λµ x ) ( 1 + 1 + θ 1 -θ ) r ′ + ( 1 + θ 1 -θ ) s ′ s.
In other words, (u k , w k ) belongs to V x . Including the hypothesis of metric regularity property once more, we find

dist ( u k , Φ -1 x (w k ) ) κdist ( w k , Φ x (u k ) ) κ ∥w k -w k-1 ∥ R n = κ ∥w -Ψ y,x (v k ) -[w -Ψ y,x (v k-1 )]∥ R n = κ ∥Ψ y,x (v k ) -Ψ y,x (v k-1 )∥ R n κ(λ ′ + λµ x ) ∥v k -v k-1 ∥ y .
On the other hand,

∥v k -v k-1 ∥ y ρ -1 1 d ( R y (v k ), R y (v k-1 ) ) = ρ -1 1 d ( R x (u k ), R x (u k-1 ) ) ρ -1 1 ρ 2 ∥u k -u k-1 ∥ x , which implies that dist ( u k , Φ -1 x (w k ) ) θ ∥u k -u k-1 ∥ x θ k ∥u 1 -u 0 ∥ x .
However, Φ -1 x (w k ) is a closed subset of a finite dimension space, so it must contain at least a vector u k+1 such that

∥u k -u k+1 ∥ x = dist ( u k , Φ -1 x (w k ) ) θ k ∥u 1 -u 0 ∥ x .
As a result,

∥u k+1 -ū∥ x k ∑ j=0 ∥u j+1 -u j ∥ x + ∥u 0 -ū∥ x k ∑ j=0 θ j ∥u 1 -u 0 ∥ x + ∥u 0 -ū∥ x 1 1 -θ ∥u 1 -u 0 ∥ x + ∥u 0 -ū∥ x .
According to [START_REF] Sakai | Riemannian Geometry[END_REF] and ( 19) it is possible to write

∥u 1 -u 0 ∥ x ∥u 0 -ū∥ x + κ ∥w 0 -p∥ R n [ρ -1 1 ρ 2 + κ(λ ′ + λµ x )] ∥v 0 -v∥ y + κ ∥w -q∥ R n < ρ -1 1 ρ 2 (1 + θ)r ′ + κs ′ . Notice that ∥u 0 -ū∥ x ρ -1 1 ρ 2 ∥v 0 -v∥ y < ρ -1 1 ρ 2 r ′ , we deduce ∥u k+1 -ū∥ x < ρ -1 1 ρ 2 ( 1 + 1 + θ 1 -θ ) r ′ + κ 1 -θ s ′ r.
In particular, u k+1 ∈ B x (ū, r), the sequences (u k ) and (v k ) are now well-defined by induction.

From the construction above, both (u k ) and (v k ) are Cauchy sequences. Thus, there are u * ∈ T x M and v * ∈ T y M such that u k → u * and v k → v * . Passing to the limit in the inclusion

w + Df (x)(u k ) -Df (y)(v k ) = w k ∈ Φ x (u k+1 ) we conclude w + Df (x)(u * ) -Df (y)(v * ) ∈ Φ x (u * ). Equivalently, w ∈ Df (y)(v * ) + (F • R x ) (u * ). But from R x (u k ) = R y (v k ) we also get R x (u * ) = R y (v * ), which gives us w ∈ Df (y)(v * ) + (F • R y ) (v * ) = Φ y (v * ). Consequently, dist ( v, Φ -1 y (w) ) ∥v -v * ∥ y = ∥v 0 -v * ∥ y = ∑ k 0 (v k -v k+1 ) y ∑ k 0 ∥v k -v k+1 ∥ y ∑ k 0 ρ -1 1 ρ 2 ∥u k -u k+1 ∥ x ρ -1 1 ρ 2 ∑ k 0 θ k ∥u 1 -u 0 ∥ x = ρ -1 1 ρ 2 1 1 -θ ∥u 1 -u 0 ∥ x ρ -1 1 ρ 2 1 1 -θ κη.
This is exactly [START_REF] Papadopoulos | Metric Spaces, Convexity and Nonpositive Curvature[END_REF]. The proof is thereby complete.

A semi-local version of Proposition 3.1 will be useful for the study of algorithm [START_REF] Adler | Newton's method on Riemannian manifolds and a geometric model for the human spine[END_REF]. The next statement is in this sense.

Proposition 3.3 (stability of semi-local metric regularity).

Let M, R, r x , W x , δ = C R (x) and ρ 1 , ρ 2 be as in the statement of Proposition 3.1. Let f : M -→ R n be a given smooth map and F : M ⇒ R n be a multivalued mapping having closed graph. Fix a point x ∈ M. Suppose that the mapping

Φ x (•) := Df (x)(•) + (F • R x ) (•) is metrically regular on the set V r,s (Φ x ) := { (u, w) ∈ T x M × R n : ∥u∥ x r, dist ( w, Φ x (u) ) s }
with a modulus κ > 0. Consider some positive numbers r ′ , s ′ , σ, λ and λ ′ such that

         θ = ρ -1 1 ρ 2 κ(λ ′ + λµ x ) < 1, ρ -1 1 σ + ρ -1 1 ρ 2 r ′ + κ 1-θ s ′ < min {δ, r} , s ′ s, ρ -1 1 ρ 2 σ + ρ -1 1 ρ 2 2 r ′ + ρ 2 κ 1-θ s ′ < r x . ( 20
)
Pick some y ∈ W x so that d(y, x) σ and the conditions (i) and (ii) of Proposition 3.1 hold. Then the multivalued mapping

Φ y (•) := Df (y)(•) + (F • R y ) (•) satisfies τ ∈ REG ( Φ y , V r ′ ,s ′ (Φ y ) ) for τ = ρ -1 1 ρ 2 κ 1-θ and V r ′ ,s ′ (Φ y ) := { (v, w) ∈ T y M × R n : ∥v∥ x r ′ , dist ( w, Φ y (v) ) s ′ } . Proof. We fix some (v, w) ∈ V r ′ ,s ′ (Φ y ) with η = dist ( w, Φ y (v) ) > 0 and look for a point v * ∈ Φ -1 y (w) such that ∥v -v * ∥ x τ η. (21)
For this goal, let us set v 0 = v, and y 0 = R y (v 0 ). Since ∥v 0 ∥ y r ′ ρ -1 1 ρ 2 r ′ < δ, the fact R ∈ URC(ρ 1 , ρ 2 , δ, W x ) can be used. This yields

d(R y (v 0 ), x) d(R y (v 0 ), y) + d(y, x) ρ 2 ∥v 0 ∥ y + σ ρ 2 r ′ + σ < ρ -1 1 ρ 2 σ + ρ -1 1 ρ 2 2 r ′ + ρ 2 κ 1 -θ s ′ < r x ,
which implies y 0 = R y (v 0 ) ∈ W x . Thus, there exists (unique) u 0 ∈ δB x with R x (u 0 ) = R y (v 0 ). By setting

w 0 = w + Df (x)(u 0 ) -Df (y)(v 0 ) ∈ R n ,
we are going to verify (u 0 , w 0 ) ∈ V r,s (Φ x ). Indeed, according to Assumption 2.2

∥u 0 ∥ x ρ -1 1 d ( R x (u 0 ), x ) = ρ -1 1 d ( R y (v 0 ), x ) ρ -1 1 ( d ( R y (v 0 ), y ) + d ( y, x )) ρ -1 1 ( ρ 2 ∥v 0 ∥ y + σ ) = ρ -1 1 ρ 2 ∥v 0 ∥ y + ρ -1 1 σ, (22)
and therefore

∥u 0 ∥ x ρ -1 1 ρ 2 r ′ + ρ -1 1 σ < ρ -1 1 σ + ρ -1 1 ρ 2 r ′ + κ 1 -θ s ′ r.
In addition, one has

dist ( w 0 , Φ x (u 0 ) ) = dist ( w 0 -Df (x)(u 0 ), (F • R x ) (u 0 ) ) = dist ( w -Df (y)(v 0 ), (F • R y ) (v 0 ) ) = dist ( w, Df (y)(v 0 ) + (F • R y ) (v 0 ) ) = dist ( w, Φ y (v) ) s ′ s.
So, the inclusion

(u 0 , w 0 ) ∈ V r,s (Φ x ) is now clear. Since κ ∈ REG ( Φ x , V r,s (Φ x ) ) , we get dist ( u 0 , Φ -1 x (w 0 ) ) κdist ( w 0 , Φ x (u 0 ) ) = κdist ( w, Φ y (v) ) = κη.
Selecting

u 1 ∈ Φ -1 x (w 0 ) that ∥u 0 -u 1 ∥ x = dist ( u 0 , Φ -1 x (w 0 )
) , we conclude

∥u 0 -u 1 ∥ x κη.
With respect to the inductive step, assume that for some k 1 we have found

u 1 , . . . , u k in T x M and v 0 , . . . , v k-1 in δB y ⊂ T y M such that • R x (u j ) = R y (v j ), j = 0, 1, . . . , k -1; • u j+1 ∈ Φ -1 x (w j ) for w j = w + Df (x)(u j ) -Df (y)(v j ) and j k -1; • ∥u j -u j+1 ∥ x θ j ∥u j -u j+1 ∥ x , j k -1.
In the tangent space T x M the triangle inequality tells us

∥u j -u 0 ∥ x j-1 ∑ i=0 ∥u i+1 -u i ∥ x j-1 ∑ i=0 θ i ∥u 1 -u 0 ∥ x = 1 -θ j 1 -θ ∥u 1 -u 0 ∥ x 1 1 -θ ∥u 1 -u 0 ∥ x , j = 1, . . . , k.
Taking into account (22), ∥u 0 ∥ x ρ -1 1 ρ 2 ∥v 0 ∥ y + ρ -1 1 σ, which yields

∥u j ∥ x ∥u j -u 0 ∥ x + ∥u 0 ∥ x 1 1 -θ κη + ρ -1 1 ρ 2 ∥v 0 ∥ y + ρ -1 1 σ ρ -1 1 σ + ρ -1 1 ρ 2 r ′ + κ 1 -θ s ′ < min {δ, r} , j = 1, . . . , k.
From Assumption 2.2, we find

d ( R x (u k ), x ) ρ 2 ∥u k ∥ x ρ -1 1 ρ 2 σ + ρ -1 1 ρ 2 2 r ′ + ρ 2 κ 1 -θ s ′ < r x .
Hence, there exists v k ∈ δB y with R y (v k ) = R x (u k ). By setting

w k = w + Df (x)(u k ) -Df (y)(v k ) ∈ R n , we claim (u k , w k ) ∈ V r,s (Φ x ).
In fact, we have known ∥u k ∥ x r. Besides that, let us now estimate the quantity dist

( w 1 , Φ x (u k ) ) as follows dist ( w k , Φ x (u k ) ) ∥w k -w k-1 ∥ R n = ∥w -Ψ y,x (v k ) -[w -Ψ y,x (v k-1 )]∥ R n = ∥Ψ y,x (v k ) -Ψ y,x (v k-1 )∥ R n (λ ′ + λµ x ) ∥v k -v k-1 ∥ y .
In the preceding estimations, we used the relations

Ψ y,x (v j ) = Df (y)(v j ) -Df (x)(u j ), j = 0, 1, . . . , k.
On the other hand, we conclude from the fact R ∈ URC(ρ 1 , ρ 2 , δ, W x ) that

∥v k -v k-1 ∥ y ρ -1 1 d ( R y (v k ), R y (v k-1 ) ) ρ -1 1 ρ 2 ∥u k -u k-1 ∥ x . Hence, dist ( w k , Φ x (u k ) ) (λ ′ + λµ x )ρ -1 1 ρ 2 ∥u k -u k-1 ∥ x (λ ′ + λµ x )ρ -1 1 ρ 2 θ k-1 ∥u 1 -u 0 ∥ x (λ ′ + λµ x )ρ -1 1 ρ 2 θ k-1 κη = θ k η θ k s ′ < s. Invoking the hypothesis κ ∈ REG ( Φ x , V r,s (Φ x ) ) , we obtain dist ( u k , Φ -1 x (w k ) ) κdist ( w k , Φ x (u k ) ) κ(λ ′ + λµ x )ρ -1 1 ρ 2 θ k-1 ∥u 1 -u 0 ∥ x = θ k ∥u 1 -u 0 ∥ x .
Since the tangent space T x M is finite dimensional, the projection u k+1 ∈

Φ -1 x (w k ) of u k onto Φ -1 x (w k ) satisfies ∥u k -u k+1 ∥ x = dist ( u k , Φ -1 x (w k ) ) θ k ∥u 1 -u 0 ∥
x , and the construction goes on. So, (u k ) and (v k ) are totally defined. The rest of the proof is similar to the one of Proposition 3.1.

Remark 3.4. In the case when the set-valued part F (x) does not depend on x (i.e., F (x) ≡ K), one has

Φ y • P χ,x,y = G χ,x,y + Φ x , (23)
for any arbitrary geodesic segment χ : [0, 1] -→ M such that χ(0) = x and χ(1) = y. Inasmuch as P χ,x,y is linearly isometric, for each subset

V ⊂ T y M × R n the fact τ ∈ REG ( Φ y , V ) is equivalent to τ ∈ REG ( Φ y •P χ,x,y , V ′ ) , in which V ′ = { (u, w) ∈ T x M × R n : (P χ,x,y u, w) ∈ V } .
It is well-known from the literature that under condition imposed on Φ x and Lipschitz modulus of G χ,x,y , the right-hand side of (23) is also metrically regular. Therefore, the suppositions related to the maps Σ χ,y,x in both Propositions 3.1 and 3.3 might be omitted. And we simply need ∥G χ,x,y ∥ λ ′ < κ -1 . Of course, the corresponding parameters as well as the region of metric regularity may be different from the old ones in Propositions 3.1 and 3.3 above.

Convergence of Newton-type algorithms

We begin with some crucial hypotheses.

Standing assumptions.

Let L 1 , L 2 : [0, +∞) -→ [0, +∞) be some nondecreasing continuous function with L 1 (0) = L 2 (0) = 0 and Ω be a nonempty open subset of M. We consider the conditions below:

(A1) Given x ∈ Ω and a geodesic χ : [0, 1] -→ Ω with χ(0) = x. If (r x , C R (x)) is a normal pair, and χ([0, 1]) ⊂ B M (x, r x ) one has Σ χ,y,x ∈ L ip L 1 (ℓ(χ)) ( C R (x)B y
) , where y = χ(1) and

Σ χ,y,x = R -1 x • R y -P χ,y,x . (A2) If Θ : [0, 1] -→ M is a geodesic in Ω joining z = Θ(0) to z ′ = Θ(1)
and

G Θ,x,y := Df (x) -Df (y) • P Θ,x,y , then ∥G Θ,x,y ∥ L 2 (ℓ(Θ)).
Under the conditions (A1) and (A2) we have the following statement.

Proposition 4.1. Given z ∈ Ω and one R-normal pair (r, C) at z. Suppose that

• the retraction segment c(t) = R z (tu) belongs to B M (z, r), • B M (z, r) is a convex neighborhood at z (cf. [7, 15]). Define Λ R (f, x, u) = f ( R x (u) ) -[f (x) + Df (x)(u)] R n , (24)
then one has

Λ R (f, z, u)) ∥u∥ z ∫ 1 0 { μL 1 ( d(z, c(t)) ) + [ L 1 ( d(z, c(t)) ) + 1 ] L 2 ( d(z, c(t)) ) } dt, (25)
where μ = ∥Df (z)∥.

Proof. Denoting Θ t , χ t : [0, 1] -→ M the minimizing geodesics with Θ t (0) = z, Θ t (1) = R z (tu), χ t (0) = z and χ t (1) = R z (tu). For simplicity we will use the notation c t = c(t), Q t = (dR z ) tu , P t = P Θt,z,c(t) , P t = P χt,z,c(t) . Then

Σ Θt,z,ct = R -1 c(t) • R z -P t and G Θt,z,ct = Df (x) -Df (c t )P t . Note that tu = R -1 z (c(t)) ∈ CB z , from (A1) and (A2) one gets d (Σ Θt,z,ct ) tu L 1 ( ℓ(Θ t ) ) , ∥G Θt,z,ct ∥ L 2 ( ℓ(Θ t ) ) .
As a consequence of the choice of

Θ t it holds that ℓ(Θ t ) = d ( z, c (t)
) . Combining this with the following relation

d (Σ Θt,z,ct ) tu = ( dR -1 c(t) ) c(t) • (dR z ) tu -P t = Q t -P t , we obtain Q t -P t L 1 ( ℓ(Θ t ) ) = L 1 ( d ( z, c(t) ))
.

Thanks to the expressions G χt,z¯,ct = Df(z¯) -Df(c t )P χt,z¯,ct , we obtain

∥Df (c t )∥ ∥-G χt,z,ct + Df (z)∥ (P χt,z,ct ) -1 ∥G χt,z,ct ∥ + ∥Df (z)∥ L 2 ( ℓ(χ t ) ) + ∥Df (z)∥ = L 2 ( d ( z, c(t) )) + μ.
The map h = f • c is smooth and

h ′ (t) = Df (c t ) ( Q t (u) ) = Df (c t ) ( Q t (u) -P t (u) ) + Df (c t ) ( P t u ) = Df (c t ) ( Q t (u) -P t (u) ) -G Θt,x,ct (u) + Df (x)(u).
Consequently,

Λ R (f, z, u) = ∥h(1) -h(0) -Df (x)(u)∥ = ∫ 1 0 h ′ (t) dt -Df (x)(u) ∫ 1 0 ( ∥Df (c t )∥ Q t -P t + ∥G Θt,x,ct ∥ ) ∥u∥ z dt.
Hence, the conclusion of Proposition 4.1 follows.

We are now ready to present the main theorems of this section. (iii) the mapping

Φ * := Df (x * )(•) + (F • R x * ) (•) satisfies τ ∈ REG ( Φ * , V * ) where V * = rB x * × B R n ( -f (x * ), s ) ; (iv) L 1 and L 2 are of class C 1 ; (v) 2K * r 1, with        K * = (ρ -1 1 ρ 2 ) 2 { µ * K 1 (r) + [L 1 (r) + 1]K 2 (r) } , K j (r) = sup 0 t r L ′ j (t) , j = 1, 2, µ * = ∥Df (x * )∥ .
We also assume in addition that W * is a convex neighborhood at x * . Set

σ = min { ρ 1 ρ 2 r, ρ 1 r, ρ 1 τ * 2K * r s, ρ 1 1 + 2K * r C * , ρ 2 2ρ 1 K * r r * } > 0, (26)
then for any x ∈ B M (x * , σ), there is a sequence (x k ) generated by scheme (2) which starts at x 0 = x and converges quadratically to x * . More precisely, one has

d(x k+1 , x * ) ρ 2 2ρ 1 r [ d(x k , x * ) ] 2 , k = 0, 1, . . . (27)
Under the suppositions of Theorem 4.2

L j (t) L j (0) + K j (r)t = K j (r)t, j = 1, 2 (28)
for every t in the interval [0, r].

Proof. Pick x ∈ B M (x * , σ) and we set x 0 = x, σ 0 = d(x * , x 0 ) < σ. If x = x * , then the proof is trivial. We focus on the case where x ̸ = x * . To find the next iteration, we shall make use of Proposition 3.1. Let χ 0 : [0, 1] -→ W * be a minimizing geodesic which links x * = χ 0 (0) and x 0 = χ 0 (1). The hypothesis (ii) of Theorem 4.2 tells us

Σ χ 0 ,x 0 ,x * ∈ L ip L 1 (σ 0 ) (W * ), ∥G χ 0 ,x * ,x 0 ∥ L 2 (σ 0 ). Choose ū = 0 x * , p = -f (x * ), v = R -1 x 0 (x * ) ∈ C * B x 0 and q = -f (x * ) + Df (x 0 )(v). Using condition (i) of Theorem 4.2 ∥v∥ x 0 ρ -1 1 d ( R x 0 (v), x 0 ) = ρ -1 1 d(x * , x 0 ) = ρ -1 1 σ 0 < ρ -1 1 σ r. Since d ( x * , R x 0 (tv) ) = d ( R x 0 (v), R x 0 (tv) ) ρ 2 ∥(1 -t)v∥ x 0 ρ 2 ρ -1 1 σ r whenever 0 t 1, the retraction segment Θ 0 (t) = R x 0 (tv) lies into W * . Moreover, d ( x 0 , Θ 0 (t)
) ρ 2 t ∥v∥ x 0 by assumption (i). Thus, applying Proposition 4.1 for z = x * and z = x 0 Λ R (f, x 0 , v)

∥v∥ x 0 ∫ 1 0 { µ * L 1 ( ρ 2 t ∥v∥ x 0 ) + [ L 1 ( ρ 2 t ∥v∥ x 0 ) + 1 ] L 2 ( ρ 2 t ∥v∥ x 0 ) } dt = 1 ρ 2 ∫ ρ 2 ∥v∥ x 0 0 { µ * L 1 (t) + [L 1 (t) + 1] L 2 (t) } dt 1 ρ 2 ∫ ρ -1 1 ρ 2 σ 0 0 { µ * L 1 (t) + [L 1 (t) + 1] L 2 (t) } dt.
Taking into account ρ -1 1 ρ 2 σ 0 < ρ -1 1 ρ 2 σ r, we get L j (t) K j (r)t for j = 1, 2

and t ∈ [0, ρ -1 1 ρ 2 σ 0 ]. As a result, Λ R (f, x 0 , v) 1 ρ 2 ∫ ρ -1 1 ρ 2 σ 0 0 { µ * L 1 (t) + [L 1 (t) + 1] L 2 (t) } dt 1 ρ 2 { µ * K 1 (r) + [L 1 (r) + 1]K 2 (r) } ∫ ρ -1 1 ρ 2 σ 0 0 t dt = ρ 2 ρ 2 1 { µ * K 1 (r) + [L 1 (r) + 1]K 2 (r) } σ 2 0 = 1 2ρ 2 τ * K * σ 2 0 .
Hence

∥f (x * ) -f (x 0 ) -Df (x 0 )(v)∥ R n = Λ R (f, x 0 , v) 1 2ρ 2 τ * K * σ 2 0 . (29)
Let us now test the following evaluations

                       θ 0 = ρ -1 1 ρ 2 τ * [L 2 (σ 0 ) + L 1 (σ 0 )µ * ] < 1, τ * 1-θ 0 [ 1 2ρ 1 τ * K * σ 2 0 ] < r, ( 1
+ θ 0 1-θ 0 ) [ 1 2ρ 1 τ * K * σ 2 0 ] < s, ρ -1 1 σ 0 + τ * 1-θ 0 [ 1 2ρ 1 τ * K * σ 2 0 ] < C * , ρ 2 τ * 1-θ 0 [ 1 2ρ 1 τ * K * σ 2 0 ] < r x . (30)
Indeed, the first relation in ( 30) is valid from the hypothesis, because of

θ 0 ρ -1 1 ρ 2 τ * σ 0 [K 2 (r) + K 1 (r)µ * ] K * σ 0 K * r 1 2 .
Besides that, one has

K * σ 0 K * r, 1 1-θ 0 1 1-K * r 2 and θ 0 1-θ 0 K * r 1-K * r 1.
So, (30) follows from the suppositions (26) and (v).

Returning to the main proof: From (30), we can take some constants r 0 > 0 and s

0 > 1 2ρ 1 τ * K * σ 2 0 with                  ρ -1 1 ρ 2 ( 1 + 1+θ 0 1-θ 0 ) r 0 + τ * 1-θ 0 s 0 < r, ( ρ -1 1 ρ 2 ) 2 [L 2 (σ 0 ) + L 2 (σ 0 )µ * ] ( 1 + 1+θ 0 1-θ 0 ) r 0 + ( 1 + θ 0 1-θ 0 ) s 0 < s, ρ -1 1 σ 0 + ρ -1 1 ρ 2 ( 1 + 1+θ 0 1-θ 0 ) r 0 + τ * 1-θ 0 s 0 < C * , ρ 2 [ ρ -1 1 ρ 2 ( 1 + 1+θ 0 1-θ 0 ) r 0 + τ * 1-θ 0 s 0 ] < r x .
Now, applying Proposition 3.1 for data ū, p, v and q above, the mapping

Φ 0 (•) = Df (x 0 )(•) + (F • R x 0 ) (•) satisfies the relation τ 0 ∈ REG ( Φ 0 , V 0 ) , where τ 0 = ρ -1 1 ρ 2 τ * 1-θ 0 and V 0 = B x 0 (v, r 0 ) × B R n (q, s 0 ). If we pick w 0 = -f (x 0 ), then w 0 -q = f (x * ) -[f (x 0 ) + Df (x 0 )(v)] .
It follows from (29) that (v, w 0 ) ∈ V 0 . Thus, we deduce

dist ( v, Φ -1 0 (w 0 ) ) τ 0 dist ( w 0 , Φ 0 (v) ) τ 0 ∥w 0 -q∥ R n = Λ R (f, x 0 , v) ρ -1 1 ρ 2 τ * 1 -θ 0 [ 1 2ρ 2 τ * K * σ 2 0 ] ρ -1 1 K * σ 2 0 .
From this, the closed set Φ -1 0 (w 0 ) contains a vector v 0 such that

∥v -v 0 ∥ x 0 = dist ( v, Φ -1 0 (w 0 ) ) ρ -1 1 K * σ 2 0
. By the triangle inequality in T x 0 M, we find

∥v 0 ∥ x 0 ∥v -v 0 ∥ x 0 + ∥v∥ x 0 ρ -1 1 K * σ 2 0 + ρ -1 1 σ 0 ρ -1 1 ( K * σ + 1 ) σ 0 < ρ -1 1 ( K * r + 1 ) σ < C * . Consequently, d(x 1 , x * ) = d ( R x 0 (v 0 ), R x 0 (v) ) ρ 2 ∥v 0 -v∥ x 0 ρ 2 ρ -1 1 K * σ 2 0 ρ 2 2ρ 1 r [d(x 0 , x * )] 2 .
Taking into account σ ρ 1 ρ 2 r and d(x 0 , x * ) < σ, we get ρ 2 2ρ 1 r d(x 0 , x * ) 1 2 . This tells us x 1 ∈ B M (x * , σ). Therefore, we can use x 1 as a new starting point, and go on.

Repeating this procedure, the sequence (x k ) generated by ( 2) is well-defined for which (27) is fulfilled. Exploiting (27) many times, we arrive at

d(x k , x * ) [ ρ 2 2ρ 1 r d(x 0 , x * ) ] 2 k -1 d(x 0 , x * ) ( 1
) 2 k -1 d(x 0 , x * ), 2
which means x k → x * quadratically. Hence, the proof is done.

In the statement of Theorem 4.2, the behavior of f and F around the solution

x * is a key point. There are many situations in practice where the informations around x * are not available. So, a global version might be quite useful and significant. The next theorem will be in this sense.

Theorem 4.3 (global analysis).

Let M, R, f and F as in Theorem 4.2.

Given a point x ∈ M and let (r x , C x ) be a R-normal pair at x so that

W x = B M (x, r x ) is a convex neighborhood. We require in addition that R ∈ URC ( ρ 1 , ρ 2 , C x , W x ) with 0 < ρ 1 ρ 2 .
Assume the following statements are fulfilled:

(i) Φ(•) := Df (x)(•) + (F • R x ) (•) satisfying τ ∈ REG ( Φ, V r,s (Φ) ) for some τ > 0, r > 0, s > 0; (ii) dist ( 0, f (x) + F (x) ) < min { s, 1 2ρ 1 τ C x } ;
(iii) the assumptions (A1) and (A2) are valid with some C 1 functions L 1 , L 2 ;

(iv) α = 2Kβ 1, where

       K = ( ρ -1 1 ρ 2 ) 2 τ { µ x K 1 (r) + [L 1 (r) + 1]K 2 (r) } > 0, β = ρ 1 τ dist ( 0, f (x) + F (x) ) , K 1 (r) = sup 0 t r |L ′ 1 (t)| , K 2 (r) = sup 0 t r |L ′ 2 (t)| , and µ x = ∥Df (x)∥; (v) ηβ min { ρ 2 1 ρ -1 r, ρ 1 ρ -1 2 r, ρ 2 1 ρ -1 2 C x , ρ 2 1 ρ -2 2 r x } , for η = 2 1+ √ 1-α . Then, there exists a solution x * of problem 0 ∈ f (x) + F (x) such that d(x, x * ) ρ 2 ρ -1 1 βη r. ( 31
)
For such a solution, it has a selection of sequence (x k ) generated by (2) and converges to x * . More precisely, one has

{ d(x k , x * ) ρ 2 ρ -1 1 ( 4 √ 1-α α λ 2 k 1-λ 2 k ) β, if α < 1, d(x k , x * ) ρ 2 ρ -1 1 2 -k+1 β if α = 1, (32)
where λ = 1-

√ 1-α 1+ √ 1-α 1.
It is clear to see that the case β(τ, x) = 0 is trivial. From now, we assume β > 0.

In order to prove Theorem 4.3, the next technical lemma will be essential.

Lemma 4.4. Let us consider the quadratic polynomial ω(t) = 1 2 Kt 2t + β. The Newton iterative sequence applied for the equation ω(t) = 0 reads t 0 = 0, t k+1 = t kω(t k ) -1 ω(t k ).

(33)

Then (t k ) is well-defined, strictly increasing and converges to the smallest root

t * = 1- √ 1-α K
of ω. Furthermore, the error bounds below are inherited:

• if α < 1, then    t * -t k 4 √ 1-α α λ 2 k 1-λ 2 k (t 1 -t 0 ) = 4 √ 1-α α λ 2 k 1-λ 2 k β, 2(t k+1 -t k ) 1+ √ 1+4λ 2 k (1+λ 2 k ) -2 t * -t k λ 2 k-1 (t k -t k-1 ); (34)
• if α = 1, then { t * -t k 2 -k+1 (t 1 -t 0 ) = 2 -k+1 β, 2 (√ 2 -1 ) (t k+1 -t k ) t * -t k t k -t k-1 . ( 35
) dist ( 0 x k , Φ -1 k (z k ) ) τ k dist ( z k , Φ k (0 x k ) ) ρ -1 1 ρ 2 τ 1 -Kt k ρ 2 2ρ 2 1 { µ x K 1 (r) + [L 1 (r) + 1]K 2 (r) } δ 2 k-1 = ρ -1 1 ( 1 1 -Kt k 1 2 Kδ 2 k-1 ) = ρ -1 1 [ -ω ′ (t k ) -1 ω(t k ) ] = ρ -1 1 (t k+1 -t k ). Because T x k M is finite dimensional, there exists v k ∈ Φ -1 k (z k ) with ∥v k ∥ x k = dist ( 0 x k , Φ -1 k (z k )
) ρ -1 1 (t k+1t k ). And we update now x k+1 = R x k (v k ) as the next iteration. The selection of v k tells us -f (x k ) = z k ∈ Φ k (v k ), which is equivalent to

-f (x k ) ∈ Df (x k )(v k ) + (F • R x k ) (v k ).
As a consequence, (x k , v k ) is generated by [START_REF] Adler | Newton's method on Riemannian manifolds and a geometric model for the human spine[END_REF]. Using (38) once more, from hypothesis (ii) we arrive at (t j+1t j ) = ρ 2 ρ -1 1 (t *t k ) .

∥v k ∥ x k 1 + √ 2 2 τ ρ 1 dist ( 0, f (x) + F (x) ) < C x ,
To finish the proof, we claim that x * solves (1). In fact, it follows from the preceding construction that So, by removing all conditions related to L 1 and putting L instead of L 2 in both Theorems 4.2 and 4.3, we obtain Kantorovich-type results for algorithm [START_REF] Adler | Newton's method on Riemannian manifolds and a geometric model for the human spine[END_REF], which might be viewed as extensions of the corresponding ones studied in [START_REF] Wang | Extended Newton's method for mappings on Riemannian manifolds with values in a cone[END_REF].

0 ∈ f (x k ) + Df (x k )(v k ) + F (x k+1 ). (44

Concluding Remarks

In this paper, we have studied a Newton-type algorithm for solving generalized equations in a finite dimensional Riemannian manifold. We provided sufficient conditions that guaranteed the existence and the quadratic convergence of the sequence generated by the extended Newton's method. This work can be considered one of the first that studied Newton's type algorithms for set-valued of the form [START_REF] Absil | Sepulchre: Optimization Algorithms on Matrix Manifolds[END_REF], where the metric regularity concept plays an important role. Due to the large number of applications in the calculus of variations or PDE's e.g., it would be interesting to consider the infinite dimensional Riemannian manifolds modeled on a Banach or a Hilbert space.

Assumption 2 . 2 .

 22 Given a retraction R : T M -→ M, an open subset Ω ⊂

Theorem 4 . 2 (

 42 local convergence). Given a complete Riemannian manifold M of dimension m and a retraction R : T M -→ M. Let f : M -→ R n be a smooth map, and F : M ⇒ R n be a closed multivalued mapping. Let x * ∈ M be a solution of problem (1), and (r * , C * ) be a R-normal pair at x * . Assume that (i) R ∈ URC ( ρ 1 , ρ 2 , C * , W * ) for some 0 < ρ 1 ρ 2 and W * = B M (x * , r * ); (ii) the statements (A1) and (A2) hold for some real functions L 1 , L 2 together with Ω = W * ;

  k+1 , x k ) = d (R x k (v k ), x k ) ρ 2 ∥v k ∥ x k ρ 2 ρ -1 1 (t k+1t k ) .The construction of sequence (x k ) is now completed by induction. Since (t k ) is convergent, and d(x k , x k+1 ) ρ 2 ρ -1 1 (t k+1t k ), the sequence (x k ) is Cauchy. So, x k must converge to some x * in M. The error bounds in (31) and (32) follows from (34), (35) and the followingd(x k , x * ) ∑ j k d(x j , x j+1 ) ρ 2 ρ -1

)

  Fixing an index k which is large enough. Ifχ k : [0, 1] -→ W x is a minimizing geodesic joining x * = χ k (0) to x k = χ k (1), then ℓ(χ k ) = d(x * , x k ). Using the hypothesis (iii) of Theorem 4.3, ∥G χ k ,x * ,x k ∥ L 2 (ℓ(χ k )) = L 2 ( d(x * , x k ) ) K 2 (r)d(x * , x k ).

Proof of Lemma 4.4. The existence of (t k ) and the error bounds (34), (35) can be found in [START_REF] Gragg | Optimal error bounds for the Newton-Kantorovich theorem[END_REF]. For the monotonicity, let us define the function φ : t ∈ (-∞, t * ) -→ tω ′ (t) -1 ω(t). By a simple computation φ ′ (t) = ω ′′ (t)ω (t) [ω ′ (t)] 2 = K ω (t) [ω ′ (t)] 2 . The polynomial ω is positive on the interval I = (-∞, t * ), so φ ′ (t) > 0. Since t k+1 = φ(t k ), we obtain the necessary conclusion by induction. Remark 4.5. From Lemma 4.4 we have

for all k. Moreover, if α = 1, error bounds in (35) yield

When α < 1, taking into account (34), we deduce

In particular,

Now, we begin to prove Theorem 4.3.

Proof of Theorem 4.3. We shall construct by induction the sequence (x k ) satisfying

At the initial step, we set

To continue, we pick

Therefore, x 1 verifies (2). Furthermore, note that η

which means that (39) is valid for k = 0.

Passing to the inductive step, given k (k

So, by (36

We are now going to estimate the value of quantity dist

But by the induction hypothesis

This suggests using Proposition 4.1. We can see that for t ∈ [0, 1]

Invoking Proposition 4.1 and taking into account

As a result,

For the goal of applying Proposition 3.3, let us check that

Indeed, since t k < t * = ηβ r (see (36)), it holds that L 1 (t k ) K 1 (r)t k and

and the first r elation i n ( 43) i s d one. F or t he s econd o ne, w e e xpand the polynomial ω with center at t k-1

Consequently,

Next, we use (38) and combine this relation with the initial condition (ii)

Finally, since

1 ηβ r x , the last inequality in (43) is valid due to the supposition (v) of Theorem 4.3. In summary, (43) is completely fulfilled.

According to (43), we can select some positive numbers r k and s k such that

, where

Thanks to (42), the inclusion (0

Thus, we obtain

Finally, by passing to the limit in (44), we get 0 ∈ f (x * ) + F (x * ). Remark 4.6. In the proofs of Theorems 4.2 and 4.3, the existence and the convergence of a Newton sequence (x k ) require the assumptions (A1) and which are imposed on covariant derivative Df as well as the retraction R. Assumption (A2) usually concerns with the so-called "Lipschitz continuity" for Df (see, e.g., [START_REF] Ferreira | Kantorovich's theorem on Newton's method in Riemannian manifolds[END_REF][START_REF] Ferreira | Local convergence of Newton's method under a majorant condition in Riemannian manifolds[END_REF]). If the higher order covariant derivatives D k f (cf. [START_REF] Dedieu | Newton's method on Riemannian manifolds: covariant alpha theory[END_REF][START_REF] Wang | Extended Newton's method for mappings on Riemannian manifolds with values in a cone[END_REF]) are involved in the hypotheses of those theorems, we obtain (A2) from the informations on D k f . Therefore, we can achieve some new versions of Theorems 4.2 and 4.3 under conditions of type Kantorovich and/or Smale. On the other hand, if we restrict our consideration on the case F (x) ≡ C, and R = exp is the exponential, then Theorems 4.2 and 4.3 can be slightly improved. Indeed, as noticed in Remark 3.4, the stability of metric regularity property can be obtained by using only (A2). Besides that, when R = exp, the essential Proposition 4.1 holds with L 1 (t) ≡ 0 and L 2 (t) is the function in (A2) since P χ,χ(0),χ(t) (χ ′ (0)) = χ ′ (t) for any arbitrary geodesic χ. Consequently, Kantorovich-type versions of Theorems 4.2 and 4.3 like in [START_REF] Wang | Extended Newton's method for mappings on Riemannian manifolds with values in a cone[END_REF] can be recovered. Now, we keep assuming F (x) ≡ C, while (A2) is replaced by the following (A3) For each x and y in Ω so that R -1

x (y) is a singleton, one has

where T R x,y is the vector transport defined by

x (y).

We refer to [1, Section 8.1] and [START_REF] Ring | Optimization methods on Riemannian manifolds and their application to shape space[END_REF] for more details about transportation of the form T R x,y . Then, by considering f along the retraction curve χ(t) = R z (tu), the estimation (25) in Proposition 4.1 simply reads Λ R (f, z, u) ∥u∥ z


Frédéric Aitken

Nelly Bonifaci

Klaus Von Haeften

Jussi Eloranta

Theoretical modeling of electron mobility in superfluid 4 He

The Orsay-Trento bosonic density functional theory model is extended to include dissipation due to the viscous response of superfluid 4 He present at finite temperatures. The viscous functional is derived from the Navier-Stokes equation by using the Madelung transformation and includes the contribution of interfacial viscous response present at the gas-liquid boundaries. This contribution was obtained by calibrating the model against the experimentally determined electron mobilities from 1.2 K to 2.1 K along the saturated vapor pressure line, where the viscous response is dominated by thermal rotons. The temperature dependence of ion mobility was calculated for several different s olvation cavity sizes and the data are rationalized in the context of roton scattering and Stokes limited mobility models. Results are compared to the experimentally observed "exotic ion" data, which provides estimates for the corresponding bubble sizes in the liquid. Possible sources of such ions are briefly discussed.

I. INTRODUCTION

An excess electron confined in superfluid 4 He is an excellent prototype system for modeling the general aspects of nonpolar solvation and the interaction of nanoscale objects with liquids. For example, the static solvation structure of electrons in liquid helium has received a considerable amount of experimental and theoretical attention. [START_REF] Borghesani | Ions and Electrons in Liquid Helium[END_REF] These studies established that electrons reside in spherical cavities ("bubbles") in the liquid with a barycenter radius of ca. 18.5 Å at saturated vapor pressure as determined by earlier density functional theory (DFT) calculations. [START_REF] Grau | [END_REF]3 In contrast, such bubbles are considerably smaller in polar liquids, such as water, due to the strong electrostatic interactions. 4 Since the bubble size in superfluid helium is determined by a delicate balance between the electron zero-point energy and the liquid surface energy, its radius is strongly dependent on the external pressure. The quantized energy levels of an electron trapped in such a spherical well have been studied experimentally (i.e., optical transitions between 1s-1p and 1s-2p states) as a function of external pressure. [5][6][7][8][9] These observations are in excellent agreement with the earlier DFT calculations, which employ several different types of electronhelium pseudopotentials to include the correlation effects into the model. [START_REF] Grau | [END_REF]3 Due to the quantum analog of the Bernoulli force, electron bubbles may also become trapped on quantized vortex lines and rings present in superfluid helium. [START_REF] Borghesani | Ions and Electrons in Liquid Helium[END_REF]10,[START_REF] Donnelly | Quantized Vortices in Helium II[END_REF] In fact, the first experimental visualization of quantized vortex lines was based on this principle. 12,13 Dynamic phenomena of solvated electrons in superfluid helium have been mostly studied experimentally where they were subject to either an external static electric field a) E-mail: Jussi.Eloranta@csun.edu (i.e., electron mobility experiments) or ultrasound excitation to explode the electron bubbles by negative pressure waves. [START_REF] Borghesani | Ions and Electrons in Liquid Helium[END_REF]14,15 The latter experiments established that electrons can be optically detected in the liquid through such expanded bubbles. 16 The electron mobility experiments, on the other hand, have demonstrated that such measurements can yield detailed information about the dissipative interaction with the surrounding superfluid. [17][18][19][20] In the low velocity regime (a few m/s), the measurements have shown that the electron mobility is reduced by collisions with the thermal excitations present in the liquid, i.e., phonons and rotons. Three different temperature regimes have been approximately identified at saturated vapor pressure: (1) phonon limited mobility when T < 0.6 K, (2) 0.6 K < T < 1.2 K where the dissipation takes place primarily by discrete roton -electron bubble collisions and phonons, and (3) 1.2 K < T < 2.17 K where the roton collision frequency is sufficiently high such that continuum models become applicable for describing the effective interaction. [START_REF] Borghesani | Ions and Electrons in Liquid Helium[END_REF] The transition temperature between cases 2 and 3 will also be discussed further in this paper as it appears slightly different than given in Ref. 1. At higher drift velocities (some tens of m/s), a critical threshold value has been observed after which the electron mobility suddenly begins to decrease. [START_REF] Borghesani | Ions and Electrons in Liquid Helium[END_REF][21][22][23][24] This transition point has been attributed to the creation of rotons or vortex rings when the electron reaches the Landau critical velocity corresponding to the temperature, pressure, and geometry used in the experiment.

Time-of-flight measurements of the solvated electrons in an external electric field have revealed that, in addition to the normal electron bubble, a sequence of additional negative ion signals were also present. Such faster than electron "exotic ion" signals have been reported over the years by many groups [25][26][27][28][29] and various explanations for their origin have been proposed in the literature: 30,31 (1) fractional electrons where the electron wavefunction is split into different solvation bubbles, i.e., a superposition state, (2) electronically excited negatively charged He -and He - 2 , 32-37 and (3) solvated impurities (e.g., O -, O - 2 H - 2 ). However, none of the previously given explanations fully account for the experimental observations and hence the problem of identifying these exotic ions remains open.

To date, theoretical modeling of the time-dependent liquid helium response in the presence of an excess electron has received relative little attention in the literature. [START_REF] Borghesani | Ions and Electrons in Liquid Helium[END_REF] For example, the early work of Rosenblit and Jortner used a hydrodynamic model to estimate the free electron localization time in liquid helium. 38 In the superfluid phase, time-dependent DFT calculations have been used to predict the variation of the energy level separations between 1s, 2s, 1p, and 2p states following an optical excitation. 3 This calculation demonstrated the feasibility of performing picosecond-times cale optical pump-probe experiments on electron bubbles to monitor their evolution. Recently, a similar DFT model has been applied to evaluate the hydrodynamic added mass of negative and positive charge carriers in superfluid helium. 39 As the current DFT models do not include the thermal excitations that are responsible for viscous dissipation, Aitken et al. have developed a phenomenological model where the hydrodynamic radius is related to the free volume and thermodynamic state equations that describe its variation with phase, pressure, and temperature. 40,41 The state equations are calibrated to the known experimental data for the particular ion and medium under consideration. This model has been successful in describing electron mobility in supercritical and liquid helium, 40 helium gas at low temperature, 41 and positive ion mobility in liquid helium 42 and supercritical helium. 43 Although the model itself is less computationally demanding than DFT and not limited to a specific system, it obscures many of the microscopic details of ion solvation and the role of dissipative processes.

In this work, we extend the time-dependent Orsay-Trento DFT (OT-DFT) model to include the viscous response of superfluid helium present at finite temperatures. This is accomplished by transforming the classical viscous response terms from the Navier-Stokes equation to DFT by using the Madelung transformation. The overall viscous response in this model is not only determined by the bulk liquid shear viscosity but also by the variation of viscous response at the inhomogeneous bubble interface. The latter contribution is obtained by calibrating the model against the experimentally known electron mobilities in superfluid helium along the saturated vapor pressure line. Finally, mobilities for several different bubble sizes in the bulk liquid are predicted as a function of temperature, which can be applied, for example, to identify unknown ionic species based on their experimentally observed mobility.

II. THEORY

A. Density functional theory model

Previously, time-independent OT-DFT has been successfully applied to model electron bubbles in superfluid helium as a function of both pressure and temperature. [START_REF] Grau | [END_REF]3,[44][45][46] In this model, the coupled non-linear Schrödinger equations corresponding to superfluid helium OT-DFT and the solvated electron are expressed as: [START_REF] Grau | [END_REF]3,46

             i ∂ψ ∂t = - 2 2m He ∆ψ + V OT [ρ] ψ + V ext [ρ e ] ψ i ∂ψ e ∂t = - 2 2m e ∆ψ e + V ext [ρ] ψ e , (1)
where m He is the helium atom mass, m e is the electron mass, ψ is the liquid order parameter with the corresponding one-particle density given as ρ = |ψ| 2 , ψ e is the electron wavefunction with ρ e = |ψ e | 2 , V OT is the non-linear OT-DFT potential (excluding the kinetic correlation and backflow terms), 46 and the external potential functional is defined as:

V ext [ρ] =  ρ(r ′ )V ps (|r -r ′ |) d 3 r ′ (2)
where V ps represents the electron-helium pseudopotential 47 and the argument ρ corresponds to either liquid helium or electron density. The previously developed finite-temperature OT-DFT includes only the static contributions of the liquid thermal response whereas the dynamic dissipative processes (e.g., viscosity) are neglected. 45 To include the dissipation due to viscous liquid into Eq. ( 1), the connection between the OT-DFT equation (first line of Eq. ( 1)) and the hydrodynamic Euler equation must first be briefly considered. By using the Madelung transformation (i.e., ψ = √ ρ exp (iS/ )), 48 the OT-DFT equation can be written as (the corresponding continuity equation not shown),

ρ m ( ∂⃗ v ∂t + ( ⃗ v • ⃗ ∇ ) ⃗ v ) = -ρ ⃗ ∇ξ (3)
where the helium mass density is defined as ρ m = m He ρ, the liquid velocity field ⃗ v = ⃗ ∇S/m He , and the effective pressure is given by:

ξ = V OT + V ext - ∆ √ ρ 2m He √ ρ . (4)
The last term on the right hand side of Eq. ( 4) is termed quantum pressure. In the Navier-Stokes equation, additional terms are included on the right hand side of Eq. ( 3) to account for viscous response. Following the Stokes hypothesis (i.e., monoatomic gas), we set the bulk viscosity to zero and write the Navier-Stokes equation as (corresponding continuity equation not shown),

ρ m ( ∂⃗ v ∂t + ( ⃗ v • ⃗ ∇ ) ⃗ v ) = -ρ ⃗ ∇ξ + ⃗ ∇ •  η ( ⃗ ∇⃗ v + ( ⃗ ∇⃗ v ) T - 2 3 ( ⃗ ∇ • ⃗ v ) ⃗ 1 ) , (5)
where η = η(r) represents the shear viscosity and ⃗ 1 is the unit tensor. The term contained between the brackets corresponds to the symmetrized linear stress tensor. Note that this equation is valid in the presence of both liquid compression ( ⃗ ∇ • ⃗ v 0) and rotational flow ( ⃗ ∇ × ⃗ v 0). To recast the viscous terms of Eq. ( 5) into the appropriate form for OT-DFT, the viscous response functional, ξ NS , must be solved from

-ρ∇ξ NS = ⃗ ∇ •  η ( ⃗ ∇⃗ v + ( ⃗ ∇⃗ v ) T - 2 3 ( ⃗ ∇ • ⃗ v ) ⃗ 1 ) . (6)
This is equivalent to the following Poisson equation for ξ NS :

∆ξ NS = -⃗ ∇ • ( 1 ρ ⃗ ∇ •  η ( ⃗ ∇⃗ v + ( ⃗ ∇⃗ v ) T - 2 3 ( ⃗ ∇ • ⃗ v ) ⃗ 1 ) ) . (7)
We assume that both ρ and ⃗ v vary slowly around the electron bubble such that the viscosity becomes a functional of density only. An often employed form for this functional dependence at a given temperature is η [ρ] ∝ ρ α where α is a constant defining the variation of viscosity across the interface ("interfacial viscosity coefficient"). While no rigorous justification for employing this form for liquids exists, it may be considered to bear similarities with the polynomial-type expansions used for expressing viscosity of gases in terms of density. 49,50 The following ansatz is employed to represent the interfacial viscosity:

η T [ρ,T] = ( ρ ρ 0 (T) ) α(T ) η 0 (T), (8)
where η 0 is the shear viscosity of the bulk liquid, which is determined by the product of the normal fluid viscosity and the normal fluid fraction (i.e., η 0 = (ρ n /ρ 0 ) η n ). Note that this power form can also approximately describe the density dependence of viscosity for classical liquids. 51,52 Since there is no known analytic form for function α(T), it must be obtained by calibration against the experimentally available electron mobility data. 14 The following empirical form was used to approximate the temperature dependence:

α(T) = a 0 + a 1 exp (a 2 T) , (9)
where the empirical constants a 0 , a 1 , and a 2 were found to be 1.73, 2.32 × 10 -10 , and 11.15 K -1 , respectively, by minimizing the difference between the viscous OT-DFT calculations and experimental mobility data. 14 To model a moving electron in superfluid helium, the first equation in Eq. ( 1) is solved in the coordinate frame moving with the electron by introducing the following kinetic energy constraint operator, Tc , to the Hamiltonian: 53 Tc

ψ = -⃗ v 0 • ⃗ pψ, ( 10
)
where ⃗ v 0 is the velocity vector corresponding to the electron movement (liquid flowing in opposite direction) and ⃗ p is the momentum vector operator for the liquid. For simplicity, the electron motion is restricted along the x-axis and hence ⃗ v 0 = (v 0 , 0, 0). Note that when periodic boundary condition is imposed on the simulation box, v 0 must be an integer multiple of 2π N x ∆x , where N x is the number of discretization points along the x-axis and ∆x is the spatial step length. Therefore, the lowest applicable constrained velocity is dictated by the spatial extent of the grid. For the present calculations, v 0 = 2.3 m/s, which is well within the experimentally observed linear mobility regime and especially much lower than the Landau critical velocity for roton/vortex creation. [START_REF] Borghesani | Ions and Electrons in Liquid Helium[END_REF] The steady-state flow solution to Eq. ( 1) can be conveniently achieved by using the imaginary time method. 46,54 Since the calculation is carried out in the moving electron frame of reference, it converges to the steady-flow solution around the electron with no net acceleration. Thus the viscous drag force and the force due to the external electric field must be equal in strength but in opposite directions, ⃗ F drag = -⃗ F field . The drag force can be evaluated from the steady-state solution of Eq. ( 1) by

⃗ F drag = -   ⃗ ∇ r ρ e (r)V ps (|r -r ′ |) ρ(r ′ )d 3 r ′ d 3 r. ( 11
)
Since the viscous response term in Eq. ( 7) introduces dissipation near the gas-liquid interface, asymmetric variations in both the liquid velocity and density appear between the front and the back of the bubble relative to the incident flow. According to Eq. ( 11), such an asymmetric density profile is responsible for producing the drag force. The incident field strength at the electron, V , is directly proportional to the external force, ⃗ F field , acting on it:

V = ⃗ F field /e
, where e is the elementary charge. In the linear regime, the electron mobility, µ, can now be calculated from (in units of m 2 V -1 s -1 ),

µ ≡ v 0 V = ev 0 ⃗ F field = ev 0 ⃗ F drag . ( 12
)
Furthermore, the following derived quantities are useful for analyzing the results: 39,[START_REF] Batchelor | Introduction to Fluid Dynamics[END_REF]

m add =  ρ(r)v x (r)d 3 r/v 0 , (13)
R = e 4π µη , (14)
where m add is the hydrodynamic added mass of the electron bubble and R is the Stokes radius of the bubble based on the Stokes law. Note that the viscosity used above corresponds to the viscosity of the whole fluid, which is determined by the product of the normal fluid fraction and the normal fluid viscosity.

To predict the dependence of mobility on the bubble radius, the following artificial external potential, V e , was employed:

V e (r) = c 0 exp (-c 1 (r -c 2 )) , (15)
where constants c 0 and c 1 where fixed to 3.8003 × 10 5 K and 1.6245 Å -1 , respectively. The bubble radius corresponding to this exponentially repulsive potential can be varied by changing the value of c 2 . For example, choosing c 2 = 10.1 Å, results in a density profile that resembles the electron bubble structure very closely. In such calculations, only the liquid portion of Eq. ( 1) was considered with V ext given by Eq. ( 15) and integration over the (non-existent) electronic degrees of freedom was omitted in Eq. (11).

The model defined by Eq. ( 1) consists of two coupled non-linear Schrödinger equations (i.e., superfluid helium and electron), which were solved numerically in imaginary time as described previously. [START_REF] Grau | [END_REF]3,46,54,[START_REF] Lehtovaara | Implementation of the Orsay-Trento density functional for describing superfluid helium[END_REF] The applied imaginary time step for the liquid degrees of freedom was 100 fs whereas a 0.1 fs time step was used for the electron due to its lower mass. The spatial grid applied in the calculations consisted of 1024 × 512 × 512 grid points with a step size of 0.8 bohr. Since the interaction between the electron and the liquid involves exchange of long wavelength phonons, the simulation box must be sufficiently large to accommodate such excitations. Note that the energy of these phonons is very low and therefore the mobility converges rather slowly in imaginary time, particularly near 2.1 K. The convergence criterion in the present calculations was based on the observed difference between the consequent mobilities calculated at every 500 imaginary time iterations. When this difference fell below 0.1% of the mobility, the calculation was considered to be converged. Typically 10 4 imaginary iterations were required to reach this threshold. However, by using a solution obtained previously at a nearby temperature or bubble size, the convergence time could be reduced significantly. Furthermore, the applied predictor-corrector scheme used for solving Eq. ( 1) allowed the imaginary time step length to be increased up to 100 fs without introducing any time step bias from the applied operator splitting scheme into the solution. 46 The bulk liquid density was maintained during the imaginary time iterations by introducing a constant chemical potential to Eq. ( 1), which corresponds to the applied OT-DFT functional. Because the calculations were carried out in imaginary time, it was possible to employ regular periodic boundary conditions in the simulation box. On the contrary, real time calculations would require the use of absorbing boundaries to eliminate back reflections of waves that originate from the accelerating object. To compute the liquid viscous response functional, the Poisson equation given by Eq. ( 7) was solved numerically in the Fourier space by using the centered finite difference formula [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] where a constant of 10 -9 a.u. was added to ρ to avoid numerical instabilities arising from the division by ρ. Superfluid helium shear viscosity values at the considered temperatures along the saturated vapor pressure curve were taken from the literature. 14 Finally, we note that the OT-DFT part of Eq. ( 1) could also be written separately for the superfluid and normal fluid order parameters ("two-fluid DFT"). In this case the OT-DFT potential should be evaluated by using the total liquid density, which is the sum of the normal and super fractions, ρ = ρ s + ρ n . Only the normal fluid equation should include the viscous response determined by η n . However, this model allows unrestricted spatial separation of the normal and superfluid components, which is not compatible with the original thermal OT-DFT functional. For the solvated electron considered in this work, both Eq. ( 1) and this two-fluid OT-DFT model produced very similar results for the calculated electron mobilities where in the latter case only a small spatial separation between the two liquid fractions was observed (<2% of bulk density). Within this model, the normal fluid accumulates in the front of the electron bubble whereas the superfluid tends to fill the void behind the bubble.

III. RESULTS AND DISCUSSION

Plotting the inverse of bulk liquid shear viscosity against the experimentally observed mobility indicates nearly linear behavior down to ca. 1.2 K as demonstrated in Fig. 1. Such behavior is predicted by the Stokes law, Eq. ( 14), assuming that the bubble radius is independent of temperature. However, careful examination shows that this result holds only very approximately and, thus, for the Stokes law to hold exactly, FIG. 1. Experimental electron mobility (µ) is plotted vs. the inverse bulk shear viscosity (η -1 0 ) as black circles. 14 The blue continuous line demonstrates the near linear relationship above 1.2 K and the red dashed line is provided as a guide to the eye to highlight the non-linear behavior below 1.2 K. A red dotted vertical line marks the approximate transition point, which separates the hydrodynamic continuum and discrete roton collision regimes (see text for discussion).

R must be temperature dependent, R = R(T). Below 1.2 K temperature the behavior becomes notably non-linear (cf. "transition point" indicated in Fig. 1), which can be interpreted as the point where the hydrodynamic continuum model breaks down and discrete roton scattering events at the electron bubble dominate the dissipative response ("roton gas").

An example of the superfluid helium density and the velocity field obtained from the imaginary time solution of Eq. ( 1) is shown in Fig. 2. Since the electron drift velocity is low (≈2 m/s), the velocity field appears nearly symmetric between the front and the back of the bubble. The presence of the viscous term given by Eq. ( 7) is responsible for producing the liquid flow pattern that moves around the bubble. Note that the behavior of the hydrodynamic variables ρ and ⃗ v may become problematic at the bubble interface and inside 8), the resulting electron mobility µ DFT from Eq. ( 12), the experimental electron mobility µ exp , 14 the calculated hydrodynamic added mass m add from Eq. ( 13) in units of m He , the Stokes radius R Stokes based on µ DFT (Eq. ( 14)), and the bubble radius at half bulk density value R ρ 0 /2 . the bubble itself where ρ ≈ 0 but ⃗ v 0. On the other hand, the order parameter representation handles such situations naturally because both its real and imaginary parts approach zero when ρ → 0 and the velocity does not need to be considered. This may be important in situations where high liquid velocities are present in low-density liquid regions (e.g., rapid bubble expansion dynamics). Ideally, for this reason, functionals that do not depend on liquid velocity or involve division by liquid density should be employed.

T (K) α µ DFT (cm 2 V -1 s -1 ) µ exp (cm 2 V -1 s -1 ) m add R Stokes (Å) R ρ 0 /2 (
Results from the viscous OT-DFT model calibration calculations are shown in Table I. At each temperature listed, minimization of the difference between the calculated (Eq. ( 12)) and the experimental mobilities with respect to parameter α(T) was carried out. The optimized values of α are plotted in Fig. 3 as a function of temperature. It can be seen that α varies smoothly from 5.20 at 2.1 K down to 1.58 at 1.4 K but begins to increase rapidly below 1.2 K. Note that this takes place near the same temperature where the mean free path of rotons becomes larger than the size of the electron bubble and, consequently, the Knudsen number (Kn) exceeds 1 (i.e., gaseous in terms of roton density). [START_REF] Borghesani | Ions and Electrons in Liquid Helium[END_REF] For parameter α this effect is very pronounced as its value essentially approaches FIG. 3. The exponent α, which defines the dependence of viscosity on liquid density, is plotted as a function of temperature. An exponential fit (see Eq. ( 9)) is shown for comparison by a red line. Note that the continuum approximation begins to break down below 1.2 K where the normal fluid Knudsen number (Kn) for rotons becomes greater than one (dotted black line). The dashed black line highlights the rapid increase below 1.2 K. infinity below 1.0 K. To extend the present model to this temperature regime, a Cunningham-Millikan-type correction factor should be included to provide an effective value for viscosity. [START_REF] Cunningham | Proc. R. Soc. A[END_REF] An approximate value for α(T) between 1.4-2.1 K can be obtained from Eq. ( 9), but for accurate calculations, the exact values provided in Table I should be employed instead. The viscosity values predicted by Eq. ( 8) are plotted in Fig. 4 as a function of the liquid density. At 1.4 K, the limited cold helium gas viscosity data [START_REF] Van Itterbeek | [END_REF] appears comparable to the values obtained from the viscous OT-DFT optimization calculations. However, the gas phase values should be somewhat lower because the cold gas can only carry thermal phonons and, at this temperature, thermal rotons also contribute to liquid viscosity. Unfortunately, such data is not available at other temperatures to allow for more exact comparison. For the electron bubble, such spatial density variations occur at the gas-liquid interface, which, in turn, contribute to the interfacial viscous drag through Eq. (7). Finally, a comparison between the experimental electron mobility data obtained along the saturated vapor pressure line and the results from the viscous OT-DFT model employing Eq. ( 8) is shown in Fig. 5. Note that this agreement would naturally be exact if the optimized FIG. 4. Dependence of liquid helium viscosity (η) on density (ρ) at selected temperatures as predicted by Eq. ( 8) with the values of α(T ) specified in Table I. The gas phase data were taken from Ref. 59 with the experimental viscosity values scaled by the normal fluid fraction to account for the twofluid nature of superfluid helium. Note the logarithmic scale for η. FIG. 5. Comparison between the experimental and viscous OT-DFT electron mobilities (µ) as a function of temperature (T ) along the saturated vapor pressure line. 14,60 The interfacial viscosity coefficient α(T ) was taken from Eq. ( 8). Note the logarithmic scale for mobility. values given in Table I were used. Despite the approximate form used for expressing α(T) with Eq. ( 8), the model still reproduces the experimental electron mobility values with a good accuracy.

The small deviations of the experimental data from the Stokes law in the linear regime (see Fig. 1) imply that the Stokes radius must include some additional temperature dependent contribution. This is demonstrated in Fig. 6 where the calculated Stokes radius reaches a maximum value of ca. 35 Å between 1.4 and 1.6 K. In general, below 2.1 K, the Stokes radius appears much larger than the ρ 0 /2 radius (≈18.5 Å) obtained from the OT-DFT density profiles. This means that the Stokes radius should be interpreted as an effective radius that incorporates additional factors such as temperature dependent changes in the interfacial viscosity. Coincidentally, when the normal liquid fraction approaches one (T ≈ 2.1 K), the Stokes and ρ 0 /2 radii become equal. Ion mobilities predicted by the developed viscous OT-DFT model for selected bubble sizes at specified temperatures are shown in Fig. 7. The calculated values are observed to follow the general power-law:

µ(R) = µ e ( R e R ) γ , (16)
where µ e is the electron mobility at the given temperature with bubble radius R e (reference system) and µ and R are the corresponding quantities for the (unknown) species.

The exponent γ = 2 when the mobility is limited by roton scattering (see, for example, Refs. 31 and 61) and γ = 1 for the continuum Stokes limit. Least squares fits of Eq. ( 16) to the calculated OT-DFT mobilities with respect to the exponent γ are indicated as continuous lines in Fig. 7. In general, the obtained exponent values vary between the two limits but appear closer to the Stokes limit within the temperature range considered. While the temperature variation of γ appears systematic, it is not clear what microscopic phenomenon is responsible for this behavior. Nevertheless, the form given by Eq. ( 16) can be used to estimate ion bubble ρ 0 /2 radii (as opposed to an effective hydrodynamic radius) based on their experimentally observed mobilities. Note, however, that such an estimate assumes that (1) ion is solvated in a bubble rather than in a snowball, 62 (2) the interfacial viscosity (see Eq. ( 8)) is independent of the bubble radius, and (3) the bubble radius is large enough that the continuum model is applicable. To avoid complications from the last assumption, experiments measuring ion mobilities in superfluid helium should be carried out close to the lambda temperature. On the contrary, however, ion mobilities can be related to the FIG. 7. Predicted ion mobilities (µ) for spherical bubbles as a function of ρ 0 /2 radius (R ρ 0 /2 ) at temperatures indicated in the legend. The bubbles were defined by the repulsive potential of Eq. ( 15) where c 2 (a.u.) was set to 6.0, 10.0, 12.0, 16.0, or 19.0 to yield bubbles with the corresponding R values of 11.1 Å, 13.4 Å, 14.6 Å, 16.9 Å, and 18.6 Å at 2.1 K. The continuous lines and exponents γ shown in the legend correspond to Eq. ( 16) with the applied (µ e , R ρ 0 /2, e ) values of (0.051, 18.6), (0.066, 18.6), (0.109, 18.5), (0.187, 18.4), (0.379, 18.4), and (0.970, 18.4) in cm 2 V -1 s -1 and Å units at temperatures 2.1, 2.0, 1.8, 1.6, 1.4, and 1.2 K, respectively.

bubble radii more accurately at low temperatures (cf. Fig. 7) where the variation of mobility vs. bubble radius is more pronounced.

An example application of the data given in Fig. 7 is to estimate the bubble radii for unknown ions, such as the exotic ions, in superfluid helium. The ion mobilities of exotic ions extrapolated to 1.2 K vary between 1.4 cm 2 V -1 s -1 (Ion #8) and 3.5 cm 2 V -1 s -1 (Ion #1). 63 Furthermore, the observed "Fast Ion" species has a very high mobility, 6.0 cm 2 V -1 s -1 . Assuming a spherical geometry for the corresponding bubble, the data shown in Fig. 7 predicts that the radius for Ion #8 is ca. 13.6 Å. Although the present calculations do not extend to radii below 11 Å, Eq. ( 16) can be used to extrapolate the OT-DFT results such that Ion #4 radius is estimated to be around 10 Å. The negative helium ions, such as He * - and He * - 2 , are expected to have solvation cavities in this size regime but their limited observed lifetime in the gas phase as well as the existence of only two such ionic species are solid arguments against assigning them to these exotic ion signals (see discussion, for example, in Ref. 31; see also Ref. 32). However, they could contribute to the experimentally observed continuous background due to ionization during their flight through the cell. 63 Our preliminary ab initio electronic structure calculations at the coupled clusters level of theory (CCSD(T)) [START_REF] Werner | , version 2012.1, a package of ab initio programs[END_REF] with the basis set specified in Ref. 65 indicate that He * - 3 species is also stable in the symmetric linear geometry (bond lengths 1.22 Å). However, at present, it is not clear if this type of negative ion series, He * - n , extends any further than n = 3. Note that typical impurities in the ion mobility experiments include oxygen from air, nitrogen from air, water, and possible residues of hydrocarbons used for pumping the sample cell. The high energy process used to produce the electrons (e.g., discharge) can also dissociate and ionize these species, producing many fragments with various electron affinities. Examples of possible ions that can form include O -, O - 2 , OH -, C -, C - 2 , H -, H - 2 , and possibly even N -and N - 2 . The neutral parents of the latter two species have, however, very low electron affinities. If the previously mentioned extrapolation procedure is applied to fast ion, 63 a radius of 4 Å is obtained. Such a small bubble could be formed by a high electron affinity atom or diatomic molecule (e.g., O -or O - 2 ). Finally, we note that the present viscous OT-DFT method is not just limited to describing electron mobility in superfluid helium, but it can also model ions with a complex solvation structure. Our preliminary calculations show that, for example, the predicted ion mobilities for Ca + and Sr + at 1.2 K are 1.06 cm 2 V -1 s -1 and 0.97 cm 2 V -1 s -1 , respectively, which can be compared with the corresponding experimental values of 1.1 cm 2 V -1 s -1 and 1.2 cm 2 V -1 s -1 . 66 Both ions have a partial snowball solvation structure (i.e., hollow core with a radius of approximately 4 Å) 67 where the first rigidly bound solvation layer follows the ion. This, in turn, reduces the ion mobility and increases its apparent size in the liquid. Note that, based on the calculated hollow core sizes and the position of the first solvation shell, 67 Ca + should have a slightly higher mobility than Sr + , which is opposite to the experimental findings. 66 The only input to these calculations is the ion-helium pair potential that was obtained from earlier electronic structure calculations. 67 Application of the developed method to describe mobilities of positive ion snowballs and the behavior of the model outside the saturated vapor pressure line will be considered in a future publication.

FIG. 2 .

 2 FIG. 2. Liquid density contours (ρ) and the associated velocity field (v) for an electron moving in superfluid helium at 2.1 K obtained by viscous OT-DFT model. The velocity field shown corresponds to (v xv 0 , v y , v z ) where the liquid background moves at 2.3 m/s.

FIG. 6 .

 6 FIG.6. Temperature dependence of electron bubble Stokes radius and radius defined by ρ 0 /2 density obtained from OT-DFT calculations. The lines connecting the discrete points are drawn as guides to the eye.

TABLE I .

 I Viscous OT-DFT optimized values for the exponent α(T ) in Eq. (

ACKNOWLEDGMENTS

Financial support from the National Science Foundation Grant Nos. CHE-0949057 and DMR-1205734 to J.E. and the Leverhulme Trust (Res. Grant F00212AH), the Royal Society (Int. Exchange Grant IE130173), and COST action CM 1405 (Molim) to K.V.H. are gratefully acknowledged.


This poster gives an overview of the IEA Wind Task for Wind Power Forecasting. The Operating Agent is Gregor Giebel of DTU, Co-Operating Agent is Will Shaw of PNNL. Collaboration in the task is

The main deliverables are an IEA Recommended Practice on the forecast solution selection and evaluation process, an up-to-date list of current projects and main project results, including datasets which can be used by researchers around the world to improve their own models, and a position paper regarding the use of probabilistic forecasts. Additionally, dissemination of relevant information in both the forecasters and the end-users community is paramount, e.g. through common workshops and webinars.

Participation is open for all organisations in member states of the IEA Annex on Wind Power, see ieawind.org for the up-to-date list.

This WP brings together global leaders in NWP models as applied to the wind industry to exchange information about future research areas. The emphasis is on improvements of the wind-related forecast performance of these models especially in typical rotor heights.

Two lists of up-to-date data are mentioned below (tall met masts and experiments). Additionally, this WP verifies and validates the improvements through a common data set to test model results upon and discuss at IEA Task meetings.

RP Forecast Selection Process NWP Improvements

In the second WP, a Recommended Practice (RP) guiding the selection process of forecasting solutions was developed: The first part, the "Forecast Solution Selection Process" deals with the selection and background information necessary to collect and evaluate when developing or renewing a forecasting solution for the power market. The second part "Benchmarks and Trials", of the series offers recommendation on how to best conduct benchmarks and trials. The third part "Forecast Evaluation", provides information and guidelines regarding effective evaluation of forecasts, forecast solutions and benchmarks and trials.

Activities

Minute Scale Forecasts Information Portal

In June 2018, IEA Wind Task 32 Lidars and 36 Forecasting held a combined workshop on Very Short Term Forecasting of Wind Power. The main tools employed were lidars, radars and SCADA data. Main results were:

• Forecasts on the minute time scale are getting more important in high-wind-penetration power systems. • A combination of weather models and instrumentation provide important information when persistence fails, namely at fast changing weather conditions, ramping and high speed wind events.

Use of Uncertainty Forecasting (WP3)

The third WP surveyed the current state of use of forecast uncertainties by the power systems sector and documents and publishes results in a report and publications.

This WP also provides outreach to users of forecasts via webinars or other means to enhance their knowledge and ability to use all available information for operations.

Phase II

The second phase of the Task adds new targets:

 Discussion of possible parts of the forecasting processes to be standardised in the future.

 Online verification and benchmarking of current NWP models with met mast data.

 Detailed review of uncertainty propagation through the modeling chain.

 Assessment of the value of probabilistic forecasts, both on a theoretical basis as well as through description of actual use cases and examples.

Figure 1 :

 1 Figure 1: The instrumentation of WFIP2, in the Northwest of the USA. Source: Joel Cline.

Figure 2 :

 2 Figure 2: Overview of a simple decision support scheme illustrating common difficulties when deciding for or against trials or common procurements. Cost, validity and output of trials are often over estimated In their usefulness, because fair evaluation requires a lot of resources, and complex problem solving can often not be verified by simple tests. A guideline for decision making is therefore under preparation by the task.



  Development of an IEA Recommended Practice for the requirements of data and instrumentation for real-time forecasting.

performance/ incentive scheme

	Mob: +45 4056 5095
	Activities	News: Recommended Practice finalised and published!
						Iinitial
				yes	forecast system plan without	no	IT
					established IT infra-	infrastructure
		Internal analysis of the requirements:	structure ?	yes	built up for multiple vendors ?	no
		Long term		Short term	Complexity level high ?	yes
		plan			plan
						no
		Run a full-scale pilot project:  involve all departments  involve external parties, stakeholders etc. Establish system requirements  Consider pilot and interim solution	     	12 months term after appox. establish basic requirements start with an interim system with 1 vendor Selection through a RFI Collect experience Set interim system and plan for long-up 18-24 months	Time constraints for selection/ contract yes no Purpose: selection of new vendor yes no k Benchmar Trial	RFP: Forecast Type & Method  Requirement list from implementation or past experience Vendor capabilities  Check experience: staff CV, references  evaluate services:  Price versus value  Support structure  Escalation structure  Redundancy structure
		IT system
		established			Establish IT system
		in pilot ?	no		in accordance to
					requirements and
		yes			experience	contract award Vendor selection /
					Results



Set up



IEA Wind Task 36 Forecasting -Phase II

	G. Giebel	W. Shaw H. Frank C. Draxl	P. Pinson	J. Zack	G. Kariniotakis	C. Möhrlen R. Bessa
	DTU Wind Energy PNNL DWD	NREL	DTU Elektro UL Renewables Mines ParisTech WEPROG INESC TEC