S Chaabi

S Rigat

F Wielonsky

A boundary value problem for conjugate conductivity equations

Keywords: Boundary value problem, Riemann-Hilbert problem, Lax pair MCS: 35J25, 35Q15, 35J15

We give explicit integral formulas for the solutions of planar conjugate conductivity equations in a circular domain of the right half-plane with conductivity σ(x, y) = x p , p ∈ Z * . The representations are obtained via a Riemann-Hilbert problem on the complex plane when p is even and on a two-sheeted Riemann surface when p is odd. They involve the Dirichlet and Neumann data on the boundary of the domain. We also show how to make the conversion from one type of conditions to the other by using the so-called global relation. The method used to derive our integral representations could be applied in any bounded simply-connected domain of the right half-plane with a smooth boundary.

Introduction and main results

We seek explicit integral expressions for real-valued solutions of the Neumann/Dirichlet problem for the following two conductivity equations, considered in a bounded simply connected domain Ω of the right half-plane H = {(x, y) ∈ R 2 , x > 0}, Div (σ∇u) = 0, Div 1 σ ∇v = 0, (

for the particular case of conductivity σ(x, y) = x p , p ∈ Z. For a general conductivity σ(x, y) > 0, x ∈ H, the two equations are conjugate in the following sense. If the function u satisfies the first equation in (1.1) then the differential form -σ∂ y udx + σ∂ x udy is closed, so by applying Poincaré theorem in the simply connected Ω, there exists a function v such that

∂ x v = -σ∂ y u, ∂ y v = σ∂ x u. (1.2)
Then, v solves the second equation in (1.1). Conversely, for such a v, the differential form σ -1 ∂ y vdx -σ -1 ∂ x vdy is closed and there exists a function u that satisfies (1.2) and the first equation in (1.1). Moreover, letting

ν = (1 -σ)/(1 + σ) ∈ R,
we have that (1.2) is equivalent to the so-called conjugate Beltrami equation,

∂f = ν∂f , (1.3)
for the complex-valued function f = u + iv in Ω. This equation differs substantially from the classical and extensively studied Beltrami equation ∂f = ν∂f . The link between the conductivity equations (1.1) and (1.3) was used in [START_REF] Astala | Calderón's inverse conductivity problem in the plane[END_REF] to investigate the Dirichlet to Neumann map. For a detailed study of the conjugate Beltrami equation (1.3) see e.g. [START_REF] Baratchart | Hardy spaces for the conjugate Beltrami equation in smooth domains of the complex plane[END_REF]. When σ = 1, the equations in (1.1) reduce to the Laplace equation and (1.2) to the Cauchy-Riemann equations. When σ differs from 1, we may speak, by analogy, of a σharmonic function u and a ν-holomorphic function f . For the conductivity σ(x, y) = x p , p ∈ N * , the equations in (1.1) rewrite as

∆u + p x ∂ x u = 0, (1.4) ∆v - p x ∂ x v = 0. (1.5)
These equations are particular case of the general elliptic equation

∆u + α x ∂ x u = 0, α ∈ R, (1.6)
which is the topic of the so-called generalized axially symmetric potential theory (GASPT), see [START_REF] Weinstein | Generalized axially symmetric potential theory[END_REF]. The name comes from the relation with the Laplacian in spaces of higher dimensions. For instance, let (r, ϕ, z) denotes the usual cylindrical coordinates in R 3 , Ω be a domain in the (x, z)-plane, and Ω ⊂ R 3 be the domain obtained by rotation of Ω about the z-axis. Then, u(r, z) is a solution of ∆u(r, z) + r -1 ∂ r u(r, z) = 0 in Ω if and only if u(r, z) is harmonic in Ω . Actually, a similar relation holds for the more general equation (1.4) when p ∈ N * . Indeed, if r = (x 2 1 + • • • + x 2 p+1 ) 1/2 and U (x 1 , . . . , x p+2 ) = u(r, x p+2 ) then one checks easily that ∆U (x 1 , . . . , x p+2 ) = ∆u(r, x p+2 ) + p r ∂ r u.

When considering a domain Ω with some simple geometry, this link with the Laplacian allows one to get explicit bases of solutions via the method of separation of variables. For instance, toroidal harmonics (i.e. Legendre functions with half-integer degrees, see e.g. [START_REF] Lebedev | Special functions and their applications[END_REF]) give a complete set of solutions to (1.6) when α = ±1, cf. [START_REF] Alladio | Analysis of MHD equilibria by toroidal multipolar expansions[END_REF][START_REF] Morse | Methods of Theoretical Physics[END_REF][START_REF] Van Milligen | Expansion of vacuum magnetic fields in toroidal harmonics[END_REF]. The fact that generalized Legendre functions can be used to solve (1.6) for complex values of α appears in [START_REF] Chaabi | Analyse complexe et problèmes de Dirichlet dans le plan : équation de Weinstein et autres conductivités non-bornées[END_REF][START_REF] Chaabi | Decomposition theorem and Riesz basis for axisymmetric potentials in the right half-plane[END_REF].

The GASPT theory has been investigated by many authors, among whom Weinstein [START_REF] Weinstein | Transonic flow and generalized axially symmetric potential theory, Symposium on theoretical compressible flow[END_REF][START_REF] Weinstein | Generalized axially symmetric potential theory[END_REF][START_REF] Weinstein | On a class of partial differential equations of even order[END_REF][START_REF] Weinstein | The method of axial symmetry in partial differential equations[END_REF], Vekua [START_REF] Vekua | New methods for solving elliptic equations[END_REF], Gilbert [START_REF] Gilbert | On generalized axially symmetric potentials[END_REF][START_REF] Gilbert | Some properties of generalized axially symmetric potentials[END_REF][START_REF] Gilbert | Bergman's integral operator method in generalized axially symmetric potential theory[END_REF][START_REF] Gilbert | Poisson's equation and generalized axially symmetric potential theory[END_REF], Henrici [START_REF] Henrici | On the domain of regularity of generalized axially symmetric potentials[END_REF][START_REF] Henrici | Vekua's theory of elliptic partial differential equations with analytic coefficients[END_REF][START_REF] Henrici | Complete systems of solutions for a class of singular elliptic partial differential equations[END_REF], Mackie [START_REF] Mackie | Contour integral solutions of a class of differential equations[END_REF], Ranger [START_REF] Ranger | Some integral transformation formulae for the Stokes-Beltrami equations[END_REF]. Equation (1.6) is related to a variety of problems in physics. In particular, when p = 1, equation (1.4) can be interpreted as the linearized version of the Ernst equation in the case of a static spacetime. In [START_REF] Lenells | Boundary-value problems for the stationary axisymmetric Einstein equations: a rotating disc[END_REF], this equation was studied in the unbounded domain consisting of the exterior of the real segment (0, ρ 0 ) in the right half-plane H and was related to the relativistic gravitational field produced by a rotating disk of matter. Equation (1.5) with p = 1 is related to the behavior of plasma in a tokamak. The goal of this device, which has a toroidal geometry, is to control location of the plasma in its chamber by applying magnetic fields on the boundary. From an assumption of axial symmetry, the problem is reduced to a plane section. Then, it follows from the Grad-Shafranov equation, a secondorder elliptic nonlinear partial differential equation, see [START_REF] Blum | Numerical simulation and optimal control in plasma physics: with applications to Tokamaks[END_REF][START_REF] Shafranov | On magnetohydrodynamical equilibrium configurations[END_REF], that the magnetic flux in the vacuum between the plasma and the circular boundary of the chamber satisfies the homogeneous equation (1.5). Note, that in this instance, the conductivity equation (1.5) takes place in an annular domain, that is a doubly connected domain.

Another interesting feature of equation (1.6) is that it satisfies recurrence and symmetry relations with respect to the coefficient α. Namely, denoting by (E α ) the equation (1.6), the following holds:

u(x, y) solves (E α ) if and only if x -1 ∂ x u(x, y) solves (E α+2 ), (1.7) u(x, y) solves (E α ) if and only if x α-1 u(x, y) solves (E 2-α ), (1.8)
see e.g. [START_REF] Weinstein | On a class of partial differential equations of even order[END_REF]. Finally, note that, since equation (1.6) is elliptic, we know, by general results, see e.g. [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF], that the Dirichlet problem with datas on the boundary ∂Ω of the domain has a unique solution, which is C ∞ in Ω.

In the present paper, we will stick, for simplicity, to the case of a domain Ω being the open disk D a = D(a, 1), lying in the right half-plane H, with a > 1 a point on the positive real axis. As briefly explained at the end of Section 3, the approach used to derive our integral representations could be adapted to a general Jordan domain Ω in H bounded by a smooth curve.

The first results of our study are the following Theorems 1.1 and 1.3, which give explicit integral representations for the solutions of equation (1.6) when α is an even and odd negative integer respectively. The method of proof follows the general scheme of the socalled unified transform method, see [START_REF] Fokas | A unified approach to boundary value problems[END_REF]. It uses a Lax pair, from which one may define a function φ(z, k) depending on a spectral parameter k, and which can be characterized as the solution of a specific singular Riemann-Hilbert problem in the k-plane. Solving this Riemann-Hilbert problem leads to the seeked integral representations.

In the sequel, for z inside the disk D a , we denote by z r the point on the circle C a = C(a, 1) which has the same imaginary part as z and lies to its right. Theorem 1.1. Let u be a solution of the equation

∆u + αx -1 ∂ x u = 0, α = -2m, m ∈ N,
in the domain D a with smooth tangential and (outer) normal derivatives u t and u n on the boundary C a . Then u admits the integral representation

u(z) = - 1 π Im ˆ(z,zr) (k -z)(k + z) m J(z, k)dk + 2Re (a r ) + u(z r ), z ∈ D a , (1.9)
where integration is on the segment (z, z r ), and the quantity a r can be explicitly computed in terms of the tangential derivatives along C a of u t and u n , up to order m -1, at z r . The function J(z, k) is given by

J(z, k) = - ˆCa W (z , k), where W (z, k) is the differential form W (z, k) = (k -z)(k + z) -m-1 ((k + z)u z (z)dz + (k -z)u z (z)dz) (1.10) = (k -z)(k + z) -m-1 ((k -iy)u t (z) + ixu n (z)) ds, (1.11)
with z = x + iy and ds the length element on C a .

Remark 1.2. In the particular case α = m = 0, equation (1.6) is Laplace equation and the solutions u are simply the functions harmonic in D a . In this case, the expression in the right-hand side of (1.9) simplifies to

1 π Im ˆ(z,zr) ˆCa W (z , k)dk + u(z r ) = 1 π Im ˆ(z,zr) ˆCa u z (z) k -z dz + u z (z) k + z dz dk + u(z r ) = 1 π Im ˆ(z,zr) -2iπu z (k)dk + u(z r ) = - ˆ(z,zr) u x (x)dx + u(z r )
which is indeed u(z). Note that in the second equality we have applied Cauchy formula to the analytic function u z (z).

Theorem 1.3. Let u be a solution of the equation

∆u + αx -1 ∂ x u = 0, α = -2m + 1, m ∈ N,
in the domain D a with smooth tangential and (outer) normal derivatives u t and u n on the boundary C a . Then u admits the integral representation

u(z) = - 1 2π Im ˆCa (k -z r )(k + z r ) m J(z, k) (k -z)(k + z) dk + u(z r ), z ∈ D a , (1.12)
where integration is on the circle C a , oriented counter-clockwise. The square root in the denominator has a branch cut along the segment (-z, z). Integration starts at z r where the square root is taken to be positive, and its determination is chosen so that it remains continuous along the path of integration. The function J(z, k) is explicitly given in terms of the tangential and normal derivatives u t and u n and their derivatives of order up to m -1 along the boundary C a . It can be written as a sum,

J(z, k) = J 0 (z r , k) + ˆ W (z , k), k ∈ C a , (1.13
(z , k) = (k -z )(k + z ),
with a branch cut along the segment (-z , z ). We choose the determination of λ(z r , k) that behaves like k (resp. -k) at infinity when Im k ≥ Im z r (resp. Im k ≤ Im z r ) and then keep a continuous determination of λ(z , k) when z moves along the path of integration.

Remark 1.4. Making use of the symmetry principles (1.7) and (1.8), one deduces easily from Theorems 1.1 and 1.3 similar integral representations for the solutions of (1.4), that is in the case of a positive integer coefficient α.

The second result concerns the correspondance between the Dirichlet and Neumann data u t and u n . Having a Lax pair is equivalent to the existence of a closed differential form. By Poincaré lemma, this leads to the vanishing of an integral on a closed contour, the so-called global relation. In the present case, it will follow from results in Section 2 that ˆCa Apparently, the reconstruction of u n from (1.14), or a similar relation, is not completely straightforward. It may be conjectured that the reconstruction of u n should be possible by an integral transform, like one of the Abel type, see e.g. [START_REF] Tovbis | Semiclassical Limit of the Scattering Transform for the Focusing Nonlinear Schrödinger Equation[END_REF]. We were not able to find such an integral transform in the present case. Our method uses rather the symmetry involved in the problem and the property of a particular linear differential equation stemming from the computation of a contour integral by Cauchy formula.

[(k -z)(k + z)] α/2-1 ((k -iy)u t (z) + ixu n (z))ds = 0, k ∈ C \ (D a ∪ D -a ). (1.
In Section 2 we recall the notion of Lax pairs and compute such pairs for the equation (1.6). We also briefly discuss the use of a Lax pair, or the related global relation, in deriving an explicit correspondance between Dirichlet and Neumann data in the simple case of the Laplace equation. The study of two specific Riemann-Hilbert problems, leading to the proofs of Theorems 1.1 and 1.3, is performed in Sections 3. The proof of Theorem 1.5 is displayed in Section 4.

Lax pairs and closed differential forms

A Lax pair for a partial differential equation P (u) = 0 is a pair of ordinary differential equations, for a function φ related to u, which are compatible precisely when P (u) = 0.

The general equation (1.6) admits such a Lax pair. Indeed, writing (1.6) with respect to the complex variables z and z, we get

u zz + α 2(z + z) (u z + u z ) = 0, (2.1)
and a possible way to find a Lax pair is by rewriting (2.1) in the form

(f (z, z)u z ) z + (g(z, z)u z ) z = 0, (2.2)
where the functions f and g have to be determined. Expanding (2.2) and comparing with (2.1), we get

f z = α 2(z + z) (f + g), g z = α 2(z + z) (f + g). (2.3)
Differentiating the first equation with respect to z, the second one with respect to z, and adding, we obtain

(f + g) zz = α 2(z + z) ((f + g) z + (f + g) z ) - α (z + z) 2 (f + g),
so that f + g satisfies the adjoint equation of (2.1) (see e.g. [START_REF] Folland | Introduction to partial differential equations[END_REF]Chapter 7] for the notion of the adjoint equation). We seek solutions in the form

(f + g)(z, z) = (z + z)A(z)B(z).
Plugging that in the previous equation leads to

(z + z)A z B z = (α/2 -1)(AB z + A z B),
and thus

-(α/2 -1)A/A z + z = k = (α/2 -1)B/B z -z,
where k may be any complex number. It is a new, additional parameter in the problem, the so-called spectral parameter. Solving for the two equations, we get, as possible solutions,

A(z) = (k -z) α/2-1 , B(z) = -(k + z) α/2-1 .
In view of (2.3), we may thus choose for f and g,

f (z, z) = (k -z) α/2 (k + z) α/2-1 , g(z, z) = -(k -z) α/2-1 (k + z) α/2 .
Hence, (2.1) is equivalent to

(k + z) α/2-1 (k -z) α/2 u z z -(k + z) α/2 (k -z) α/2-1 u z z = 0.
This last equation is equivalent to the compatibility of the two ordinary differential equations

φ z (z, k) = (k + z) α/2 (k -z) α/2-1 u z (z), φ z (z, k) = (k + z) α/2-1 (k -z) α/2 u z (z), (2.4)
which thus gives a Lax pair for (2.1). Equivalently, we may express the property of the Lax pair as the closedness of a differential form.

z → W (z, k) = (k -z)(k + z) α/2-1 (k + z)u z (z)dz + (k -z)u z (z)dz (2.5)
is closed in Ω. Note that, when α ∈ 2N * , the differential form has no singularities in Ω and k can be any complex number. Otherwise, for α ∈ R \ 2N * , W (z, k) has either a pole or a branch point at k and -k, so that k should lie outside Ω and -Ω.

Making use of the symmetry relation (1.8), we derive a second Lax pair, namely

φ z (z, k) = (k + z) 1-α/2 (k -z) -α/2 ((z + z) α-1 u) z (z) = (k + z) 1-α/2 (k -z) -α/2 (z + z) α-2 [(z + z)u z (z) + (α -1)u(z)], φ z (z, k) = (k + z) -α/2 (k -z) 1-α/2 ((z + z) α-1 u) z (z) = (k + z) -α/2 (k -z) 1-α/2 (z + z) α-2 [(z + z)u z (z) + (α -1)u(z)],
and the related closed form

z → W (z, k) = (k -z)(k + z) -α/2 x α-2 • (k + z)(2xu z (z) + (α -1)u(z))dz + (k -z)(2xu z (z) + (α -1)u(z))dz . (2.6)
This second differential form will be useful in Section 4.

We end this section with a brief discussion on the possible use of a Lax pair for the derivation of an explicit correspondance between different types of boundary data. It is easy to check that not all Lax pairs can achieve this goal. For instance, if we consider (2.4) when m = 0, we get

φ z (z, k) = u z (z)/(k -z), φ z (z, k) = u z (z)/(k + z),
which is a Lax pair for the Laplace equation ∆u = 0. If we choose as a domain Ω the unit disk D, the differential form W (z, k) from Proposition 2.1 is closed in D and thus we get, in terms of the tangent and normal derivatives u t and u n on the unit circle,

ˆ∂D xu n (z) (k -z)(k + z) ds = ˆ∂D (y + ik)u t (z) (k -z)(k + z) ds, k ∈ C \ D. (2.7)
Assume that the tangent derivative u t is known, hence also the integral in the right-hand side, equal to some function ψ(k). Decomposing the real function xu n (z) = g(z) + g(1/z), with g(z) analytic in D, g(0) ∈ R, it is readily checked that (2.7) rewrites as

g(1/k) -g(-1/k) = kψ(k),
which allows one to recover only the half of g, namely the imaginary parts of its even Taylor coefficients and the real parts of the odd ones. On the contrary, when considering the domain D a = D(a, 1) instead of D, we will see in Section 4 that the Lax pair in (2.4) is sufficient to recover the Dirichlet to Neumann correspondance.

Characterization by a Riemann-Hilbert problem

The goal of this section is to prove Theorems 1.1 and 1.3.

Case of a negative even integer coefficient

α = -2m, m ∈ N Proof of Theorem 1.1. The differential form W (z, k) rewrites as W (z, k) = (k -z)(k + z) -m-1 (k + z)u z (z)dz + (k -z)u z (z)dz . (3.1)
Equations ( 2.4) can be written as dφ = W and we thus construct a function φ of the form

φ(z, k) = ˆz zr W (z , k), (3.2)
where the path of integration needs to be defined. For k such that Im k ≥ Im z, we integrate from z r to z following the lower part C low of the circle C a and then the segment from z l to z. For k such that Im k ≤ Im z, we integrate from z r to z following the upper part of the circle C up and then again the segment from z l to z, see Figure 1. This defines φ(z, k) as an In the sequel, we use the notations φ + (z, k) and φ -(z, k) for the limit values of φ(z, k) when k tends to the left and right of a given arc, according to its orientation. The fact that the differential form W (z , k) is closed in the disk D a when k lies outside of that disk implies that there is no jump

a z z l -z -a C up 0 z r C low ∆ z
J(z, k) := φ + (z, k) -φ -(z, k),
of φ(z, k) on the part of ∆ z outside of D a and -D a . For k inside the disks, we have

J(z, k) = - ˆCa W (z , k), k ∈ (z, z r ) ∪ (-z r , -z), J(z, k) = 0, k ∈ (z l , z) ∪ (-z, -z l ).
For computing the jump on (-z r , -z), we have used the symmetry (3.3). For the jump on (z l , z), we have deformed the two paths of integration into the segment (z r , z) which is possible by closedness of the differential form W (z , k). Note that k → J(z, k) has no singularity in z and z r . In z it is clear, while in z r , it can be seen by performing m integrations by parts of W (z , k) on the closed contour C a , see also [START_REF] Gakhov | Boundary Value Problems[END_REF]Section 4.4] for a similar computation. Concerning the function k → φ(z, k), as said before, it has poles of order m at {z, z r , -z, -z r } if m > 0 and logarithmic singularities if m = 0. Moreover, at infinity, it behaves like

φ(z, k) → k -2m-1 (u(z) -u(z r )), as k → ∞.
If m > 0, we renormalize the problem at infinity by defining

φ(z, k) = ((k -z)(k + z)) m φ(z, k), (3.4)
J(z, k) = ((k -z)(k + z)) m J(z, k). (3.5)
Then, the function k → φ(z, k) behaves like k -1 (u(z) -u(z r )) at infinity and is regular at z and -z, except for logarithmic singularities when m = 0. It still has polar singularities at z r and -z r . Let us denote by φ zr,-zr (z, k) its polar part at these points, that is the sum of the terms of negative degrees in the Laurent expansions of φ(z, k) at z r and -z r . The function φ -φ zr,-zr is analytic outside of the segments (z, z r ) and (-z r , -z) where it has the jump J(z, k). It has at most logarithmic singularities at {z, z r , -z, -z r } and vanishes at infinity. These properties completely determine the function φ -φ zr,-zr . Indeed, if there were another function with these properties, their difference would be entire and vanishing at infinity, hence the zero function. Thanks to the Plemelj formula, we have an integral expression for φ(z, k) -φ zr,-zr (z, k), namely,

φ(z, k) -φ zr,-zr (z, k) = 1 2iπ ˆ(-zr,-z)∪(z,zr) J(z, k ) k -k dk . (3.6)
Denoting by a r and a -r the residues of φ(z, k) at k = z r and k = -z r respectively, and equating the coefficients of k -1 in the expansion of (3.6) at infinity, we get

u(z) -u(z r ) = a r + a -r - 1 2iπ ˆ(-zr,-z)∪(z,zr) J(z, k )dk = 2Re (a r ) - 1 π Im ˆ(z,zr) J(z, k )dk ,
where we have used the symmetry relation φ(z, -k) = -φ(z, k). It remains to show that the residue a r can be explicitly computed from the knowledge of the derivatives u t and u n on C a . From (3.4), it is sufficient to know the polar part of φ(z, k) at z r . To compute this polar part, we first rewrites W (z, k), defined in (3.1), by expressing the complex derivatives in terms of the tangential and (outer) normal derivatives on the circle C a ,

u z dz = 1 2 (u t + iu n )ds, u z dz = 1 2 (u t -iu n )ds, (3.7)
with ds the length element on C a . We get, with z = x + iy ∈ C a ,

W (z, k) = (k -z)(k + z) -m-1 ((k -iy)u t (z) + ixu n (z)) ds, = (k -z) -m-1 w(z, k)dz where w(z, k) := (k + z) -m-1 ((k -iy)u t (z) + ixu n (z)) τ -1 (z)
and τ (z) denotes the unit vector tangent to C a at the point z. Let us set

∂ t f := τ -1 (z)∂ t f (3.8)
for a function f on C a . Because of analiticity,

∂ z (k -z) -j-1 dz = ∂ t (k -z) -j-1 ds = ∂ t (k -z) -j-1 dz, j = m -1, . . . , 0.
Hence, performing m integrations by parts on the integral in (3.2), we obtain

φ(z, k) = c 0 (k -z ) -m w(z , k) z zr + • • • + c m-1 (k -z ) -1 ∂ (m-1) t w(z , k) z zr -c m-1 ˆz zr (k -z ) -1 ∂ (m) t w(z , k)dz , (3.9)
where c j = (-1) j Γ(m -j)/Γ(m + 1), j = 0, . . . , m -1.

The polar part of φ(z, k) at z r can be read in the bracketed terms in (3.9) as

-c 0 m-1 j=0 ∂ (j) k w(z r , z r ) j!(k -z r ) m-j -• • • -c m-1 ∂ (m-1) t w(z r , z r ) k -z r ,
where ∂ k denotes the operator of differentiation with respect to the variable k. This finishes the proof of Theorem 1.1.

Case of a negative odd integer coefficient

α = -2m + 1, m ∈ N Proof of Theorem 1.3.
The differential form (2.5) now rewrites as

W (z, k) = (k -z)(k + z) -m-1/2 (k + z)u z (z)dz + (k -z)u z (z)dz . (3.10)
The novelty with respect to the even coefficient case is that

W (z, k) involves the square root λ(z, k) = (k -z)(k + z)
which, as a function of k, is defined on a Riemann surface S z of genus 0, consisting of two copies of C, denoted by S z,1 for the upper sheet, and by S z,2 for the lower sheet. The two sheets are glued together along a branch cut from z to -z, that we choose to be the horizontal segment (-z, z). We denote by λ 1 and λ 2 the determination of the square root λ, as a function of k, on the upper and lower sheets of S z where we assume that

λ 1 (z, k) = k(1 + O(1/k)), as k → ∞ 1 on the upper sheet S z,1 , (3.11) λ 2 (z, k) = -k(1 + O(1/k)), as k → ∞ 2 on the lower sheet S z,2 .
(3.12)

We also denote by W 1 and W 2 the values of the differential form W corresponding to the determinations λ 1 and λ 2 of the square root. For future use, we remark that the function λ(z, k) and the differential form W (z, k) satisfy the following symmetry relations:

λ(z, -k) = -λ(z, k), W (z, -k) = W (z, k). (3.13)
As in the previous section, we construct a function φ of the form

φ(z, k) = ˆz zr W (z , k), (3.14)
where the path of integration needs to be defined. Similarly to the method applied in [START_REF] Lenells | Boundary-value problems for the stationary axisymmetric Einstein equations: a rotating disc[END_REF], we define, for each z, the function k → φ(z, k) as a map from the Riemann surface S z to C. The path of integration in (3. symmetric circles, the branch cut from z to -z never intersects k as z goes from z r to z along γ 1 or γ 2 . For such k, we may therefore use the same determination of the square root λ(z , k) for computing the integral (3.14). When k lies inside one of the circles (i.e. regions I and II in Figures 2 and3), the branch cut intersects k once. Hence, for k ∈ S z,1 lying in the upper part of one of the two circles, we use the determination of the square root λ 2 (z , k) on the first part of integration, before the crossing, and the determination λ 1 (z , k) after the crossing, and conversely for k ∈ S z,2 lying in the lower part of one of the two circles. Finally, when k lies between the two circles, the branch cut intersects k twice, so that, in this case, we start integration with the determination of the square root λ(z , k) corresponding to k, then change to the other determination after the first crossing, and come back to the first determination after the second crossing. The function φ(z, k) is an analytic function of k ∈ S z outside of arcs where it has jumps. Note that, in view of (3.10), it has poles of order m at k ∈ {z, z r , -z, -z r }.

As in the previous section, φ + (z, k) and φ -(z, k) denote the limit values of φ(z, k) when k tends to the left and right of a given contour, according to its orientation. Let us compute the jumps,

J(z, k) = φ + (z, k) -φ -(z, k),
of the function φ(z, k) on the different arcs shown in Figures 2 and3. As in the previous section, we denote by C up the part of the circle C a = C(a, 1) above the segment (z l , z r ) and by C low the part of the circle C a below that segment. For k on one of the two circles, we denote by k the other point on the circle with the same imaginary part, that is, we have

k -a = -(k -a), k ∈ C a k + a = -(k + a), k ∈ -C a .
To compute the jumps J(z, k) we first rewrites W (z, k), defined in (3.10), in terms of ds, the length element on C a . By the same computation as the one in the previous section, we now get, with z = x + iy ∈ C a ,

W (z, k) = (k -z)(k + z) -m-1/2 ((k -iy)u t (z) + ixu n (z)) ds, = (k -z) -m-1/2 w(z, k)dz where w(z, k) := (k + z) -m-1/2 ((k -iy)u t (z) + ixu n (z)) τ -1 (z)
and τ (z) still denotes the unit vector tangent to C a at the point z. Performing m integrations by parts on the integral in (3.14), we obtain in a similar way as in the previous section,

φ(z, k) = c 0 (k -z ) -m+1/2 w(z , k) z zr + • • • + c m-1 (k -z ) -1/2 ∂ (m-1) t w(z , k) z zr -c m-1 ˆz zr (k -z ) -1/2 ∂ (m) t w(z , k)dz , (3.15)
where

c j = (-1) j Γ(m -1/2 -j)/Γ(m + 1/2), j = 0, . . . , m -1,
and the operator ∂ t f is still defined by (3.8). The bracketed terms contain the polar parts of φ(z, k), of degree m, at z and z r . Note that the last integral converges when k ∈ C a . Assume k ∈ S z,1 and lies on the right half of C up . Then, from the above definition of φ, we derive that

J(z, k) = J 0 1 (z r , k) + ˆk zr W 1 (z , k), with J 0 1 (z r , k) = 2 m-1 j=0 c j (k -z r ) j-m+1/2 ∂ (j) t w 1 (z r , k), (3.16
)

W 1 (z , k) = 2c m-1 (k -z ) -1/2 ∂ (m) t w 1 (z , k)dz , (3.17)
where the subscript 1 in the above expressions means that we use the determination λ 1 of the square root to evaluate them. The jumps on the left half of C up on S z,1 and on C low on S z,2 can be computed in the same way. The jumps on -C up and -C low can be derived from the symmetry relation satisfied by φ(z, k),

φ(z, -k) = φ(z, k).
The result is as follows. For k ∈ S z,1 ,

J(z, k) = J 0 1 (z r , k) + ˆ(zr,k) W 1 (z , k), k on the right half of C up , (3.18) J(z, k) = J 0 1 (z r , k) + ˆ(zr, k) W 1 (z , k) - ˆ( k,k) W 1 (z , k), k on the left half of C up , (3.19) J(z, k) = J(z, -k), k on -C up , (3.20)
where the path of integration (a, b) in each of the above integrals is the subarc of C up from a to b, in the positive direction. For k ∈ S z,2 , we have

J(z, k) = J 0 2 (z r , k) + ˆ(zr,k) W 2 (z , k), k on the right half of C low , (3.21) J(z, k) = J 0 2 (z r , k) + ˆ(zr, k) W 2 (z , k) - ˆ( k,k) W 2 (z , k), k on the left half of C low , (3.22) J(z, k) = J(z, -k), k on -C low , (3.23)
where (a, b) in each of the integrals denotes the subarc of C low from a to b, in the negative direction. Note that the above jumps take place on two open contours on the Riemann surface S z . The first one starts at z r,1 on the first sheet, follows the circle C a , goes through the cut at z l and finishes at z r,2 on the second sheet. The second one start at -z r,1 on the first sheet, follows the circle -C a , goes through the cut at -z l and finishes at -z r,2 on the second sheet. Note also that the jump on the first contour is continuous at z l (compare (3.19) and (3.22) at k = z l ), and, by symmetry, the same is true of the jump on the second contour at -z l .

It follows, simply from the definition (3.14) of φ(z, k) and the choice of the determinations, that, on each sheet, φ(z, k) has no jumps on (z, z r ) and (-z r , -z). On the branch cut (-z, z) of the Riemann surface S z there is also no jumps. Indeed, on the part (-z l , z l ), k is not close from the paths of integration and we can use the definition (3.14) to compute the jumps. One has

φ + 1 (z, k) -φ - 2 (z, k) = ˆz zr W + 1 (z , k) - ˆz zr W - 2 (z , k) = 0, φ - 1 (z, k) -φ + 2 (z, k) = ˆz zr W - 1 (z , k) - ˆz zr W + 2 (z , k) = 0,
where we have deformed the original paths of integration γ 1 and γ 2 into the segment (z r , z), which is possible since the differential form W (z, k) is closed. When k ∈ S z,1 is on theside of (z l , z) or k ∈ S z,2 is on the + side of (z l , z), we may deform both paths γ 1 and γ 2 into the segment (z r , z) and we get

φ - 1 (z, k) -φ + 2 (z, k) = ˆz zr W - 1 (z , k) - ˆz zr W + 2 (z , k) = 0.
Finally, when k ∈ S z,1 is on the + side of (z l , z) or k ∈ S z,2 is on the -side of (z l , z), we deform both paths γ 1 and γ 2 into the segment (z r , k) followed by the segment (k, z). Here k is closed from the paths of integration, so we use (3.15) to compute the jumps. We get

φ + 1 (z, k) = - 1 2 ˆk zr W - 1 (z , k) - 1 2 ˆz k W + 1 (z , k) = - 1 2 ˆk zr W + 2 (z , k) - 1 2 ˆz k W - 2 (z , k) = φ - 2 (z, k),
where we note that the bracketed terms in (3.15) do not contribute to the jumps. For completeness, let us remark that, if we would consider φ(z, k) as a function on one of the two sheets only, e.g. S z,1 , then it would have a jump on the segment (-z l , z l ), considered as an arc on S z,1 , namely ˆCup W 1 (z , k), k ∈ (-z l , z l ).

On the second sheet S z,2 , a jump would also occurs,

- ˆClow W 2 (z , k), k ∈ (-z l , z l ),
which is actually opposite to the previous one. The term J 0 (z r , k) has a polar singularity of order m at k = z r so that the jumps (3.18)-(3.23) have polar singularities of order m at either z r or -z r . Hence, instead of φ(z, k), we consider

φ(z, k) = (k -z r )(k + z r ) m φ(z, k), whose jumps J(z, k) = (k -z r )(k + z r ) m J(z, k) (3.24)
are regular (and lie on the same contours as those of φ(z, k)). From the definitions of W and φ, and in view of (3.11)-(3.12), we have lim

k→∞ 1 φ(z, k) = ˆz zr du = u(z) -u(z r ), (3.25)
lim

k→∞ 2 φ(z, k) = - ˆz zr du = u(z r ) -u(z), (3.26)
so that, in particular,

φ(z, ∞ 1 ) = -φ(z, ∞ 2 ). (3.27)
Next, the function k → φ(z, k) remains bounded near the four endpoints z r,1 , z r,2 , -z r,1 , -z r,2 of the two contours where the jumps (3.18)-(3.23) take place. Indeed, near z r,1 , W (z, k) is of order (k -z r,1 ) -m-1/2 and consequently φ(z, k) is of order (k -z r,1 ) 1/2 . The same fact holds true near the three other points. As a last remark, let us mention that φ(z, k) has poles of order m at z and -z since this holds true for φ(z, k). We denote by φ z,-z (z, k) the sum of its polar parts at z and -z.

The jumps (3.24), the relation (3.27) between the values at infinities and the boundedness near the endpoints completely characterize the function φ -φ z,-z on S z . Indeed, if there are two such functions, then their difference would be analytic on the compact Riemann surface S z , hence constant. Since it would also satisfy the relation (3.27), it could only be the zero function. An explicit expression can be given for the unique solution of the Riemann-Hilbert problem defined by the previous conditions (jumps, relation between the values at ∞ 1 and ∞ 2 , and boundedness near the endpoints), namely

φ(z, k) -φ z,-z (z, k) = 1 4iπ ˆCup∪-Cup J(z, k ) λ(z, k) λ 1 (z, k ) + 1 dk k -k + 1 4iπ ˆClow ∪-C low J(z, k ) λ(z, k) λ 2 (z, k ) + 1 dk k -k , (3.28)
where the contours of integration are oriented as in Figures 2 and3. In the first integral k ∈ S z,1 and in the second integral k ∈ S z,2 . Let us check that the expression in (3.28) satisfies the characterizing properties of φ -φ z,-z . Indeed, it defines an analytic function of k on S z outside of the two contours C up ∪ C low and -C up ∪ -C low . Thanks to the Plemelj formula, we see that it has the right jumps on these contours. Because of (3.11)-(3.12), it also satisfies relation (3.27). From the formulas (3.18), (3.20), (3.21), (3.23), we get that J(z, k ) vanishes at the endpoints of the two contours. Hence, the expression in (3.28) remains bounded near these endpoints, see [1, Lemma 7.2.2] for details. Making use of (3.25), we deduce from (3.28), where we rename k as k, that

u(z) -u(z r ) = - 1 4iπ ˆCup∪-Cup J(z, k) λ 1 (z, k) dk - 1 4iπ ˆClow ∪-C low J(z, k) λ 2 (z, k) dk = - 1 2π Im ˆCup J(z, k) λ 1 (z, k) dk - 1 2π Im ˆClow J(z, k) λ 2 (z, k) dk = - 1 2π Im ˆCa J(z, k) (k -z)(k + z) dk.
Note that the polar part φ z,-z does not give any contribution in the above computation as k tends to infinity. In the second equality we have used the first identity in (3.13) and the fact that J(z, -k) = J(z, k).

In the last expression, we integrate along the circle C a oriented counter-clockwise, starting at z r with the determination of the square root that behaves like k at infinity. At k = z l , the determination changes so that the square root remains continuous along the path of integration. This finishes the proof of Theorem 1.3.

The method used in the proofs of Theorems 1.1 and 1.3 could be applied in case of a general, bounded, simply-connected domain Ω in the right half-plane H. For the definition of the corresponding function φ(z, k), in particular the paths of integration γ 1 and γ 2 , one could replace the points z r (resp. z l ) by e.g. a point ω r (resp. ω l ) on the boundary of Ω of smallest (resp. largest) abscissa (hence independant of z). The segment (z l , z) may be replaced with any path from ω l to z. In case of an even coefficient α, another path from ω r to z can be fixed to obtain a path from ω r to ω l that separates Ω into a lower and an upper part.

The Dirichlet-Neumann map for the case of an even coefficient α

This section is devoted to the proof of Theorem 1.5. Theorems 1.1 and 1.3 give an expression of the solution to (1.5), and in particular of the jump J(z, k), in terms of its tangential and normal derivatives on the boundary of the disk D a . In practice, only one type of boundary data is usually known and it is thus important to determine if this information is sufficient for computing the solution.

Using the link (3.7) between complex and directional derivatives, the global relation, derived from Proposition 2.1 and Poincaré lemma, ˆCa

[(k -z)(k + z)] α/2-1 ((k + z)u z (z)dz + (k -z)u z (z)dz) = 0, can be rewritten as ˆCa [(k -z)(k + z)] α/2-1 ((y + ik)u t (z) -xu n (z))ds = 0, (4.1)
for k ∈ C \ (D a ∪ D -a ). We first show that, when α is an even negative integer, (4. from the relation

ˆT z m-1 f (z) (z -(k -a)) m (z + (k + a) -1 ) m dz = ˆT z m-1 (y + ik)u t (z + a) (z -(k -a)) m (z + (k + a) -1 ) m dz, (4.3)
where k ∈ C \ (D a ∪ D -a ) and we integrate on T instead of C a . Now, we set

µ := -1/(k + a), ϕ(z) := -z 1 + 2az
.

The map ϕ is involutive and sends µ to 1/(k -a). Moreover, k ∈ C \ (D a ∪ D -a ) is equivalent to the fact that µ (or ϕ(µ)) belongs to the annulus

A := D \ D -2a 4a 2 -1 , 1 4a 2 -1 , so that (4.3) can be rewritten ˆT z m-1 f (z) (1 -ϕ(µ)z) m (z -µ) m dz = 1 µ ˆT z m-1 (µy -i(aµ + 1))u t (z + a) (1 -ϕ(µ)z) m (z -µ) m dz, µ ∈ A. (4.4)
The functions f (z), u t (z+a), yu t (z+a) are real-valued in L 2 (T), so they can be decomposed as

f (z) = g 1 (z) + g 1 (1/z), u t (z + a) = g 2 (z) + g 2 (1/z), yu t (z + a) = g 3 (z) + g 3 (1/z), (4.5)
with g i (z) holomorphic on D and g i (1/z) holomorphic on C \ D for i = 1, 2, 3. Note that the imaginary parts of g 1 (0), g 2 (0), g 3 (0) can be chosen at will in the decompositions (4.5). For i = 1, 2, 3, let h i (z) := z m-1 g i (z), and

Φ i (µ) := ˆT z m-1 g i (z) (1 -ϕ(µ)z) m (z -µ) m dz = 2iπ (m -1)! h i (z) (1 -ϕ(µ)z) m (m-1) (µ), µ ∈ A.
(4.6) It is clear that the functions Φ i are analytic in A. Since 1 -ϕ(µ)µ = (µ 2 + 2aµ + 1)/(1 + 2aµ), one derives from the second expressions in (4.6) that the Φ i extend analytically to D, except at the point z 1 , where we denote by z 1 and z 2 the two roots of z 2 + 2az + 1 = 0,

z 1 = -a + √ a 2 -1 ∈ D, z 2 = ϕ(z 1 ) = -a - √ a 2 -1 ∈ C \ D.
At the point z 1 , the functions Φ i have a polar singularity of order at most 2m -1. We also note that the Φ i have a zero of order at least m at -1/2a. Next, from the fact that ϕ is involutive, follows that, for any function g(z) on T,

ˆT z m-1 g(1/z) (1 -ϕ(µ)z) m (z -µ) m dz = - ˆT z m-1 g(z) (1 -µz) m (z -ϕ(µ)) m dz.
Therefore, the relation (4.4) can be rewritten as

Φ 1 (µ) -Φ 1 (ϕ(µ)) = Φ 3 (µ) -Φ 3 (ϕ(µ)) -i(a + 1/µ)(Φ 2 (µ) -Φ 2 (ϕ(µ))), µ ∈ A. (4.7)
Note that g 2 (0) = 0 since the Fourier coefficient of u t (z + a) of order 0 vanishes. In view of (4.6), this entails that Φ 2 (0) = 0 so that the function Φ 2 (µ)/µ is analytic at 0. Multiplyimg both sides of the equation by

S(µ) = (µ -z 1 ) 2m-1 (µ -z 2 ) 2m-1 /(2aµ + 1) m ,
and rearranging terms, we get

S(µ)(Φ 1 (µ) -Φ 3 (µ) + i(a + 1/µ)Φ 2 (µ)) = S(µ)(Φ 1 (ϕ(µ)) -Φ 3 (ϕ(µ)) + i(a + 1/µ)Φ 2 (ϕ(µ))).
The function on the left-hand side is analytic in D. Moreover, the function ϕ maps A onto itself and sends the inner disk D onto C \ D with the inner circle of A onto the unit circle T. Also, -1/2a and ∞ map each others by ϕ. Hence, the term on the right-hand side is analytic in C \ D and has a polar singularity of order 2m -2 at infinity. Equating the left-hand side to the part in the right-hand side that is analytic in D, we get

S(µ)Φ 1 (µ) = S(µ)Ψ(µ) + P 2m-2 (µ), Ψ(µ) := Φ 3 (µ) -i(a + 1/µ)Φ 2 (µ),
where Ψ(µ) is a known function depending on u t , and P 2m-2 (µ) is a polynomial of degree at most 2m -2, that depends on both u t and u n , (hence is unknown). Making use of the second expressions in (4.6) of the Φ i functions and applying Leibniz rule on the (m -1)-th derivatives, we readily obtain that h 1 (z) satisfies the differential equation of order m -1,

m-1 p=0 α p z m-1-p (z 2 + 2az + 1) p h (p) 1 (z) = H(z) + P 2m-2 (z), (4.8)
where

α p = (-1) m-1+p (2m -p -2)! (m -p -1)!p! , p = 0, . . . , m -1,
and H(z) is a known function, analytic in D, that is given in terms of h 2 (z) and h 3 (z) by

H(z) = m-1 p=0 α p z m-1-p (z 2 + 2az + 1) p (h (p) 3 (z) -i(a + 1/z)h (p) 2 (z)).
By Lemma 4.1 which, for convenience, we have stated and proved after the present proof, there exists a unique solution P 2m-2 (z), polynomial of degree at most 2m -2, to the differential equation

m-1 p=0 α p z m-1-p (z 2 + 2az + 1) p P (p) 2m-2 (z) = P 2m-2 (z).
Hence, h 1 (z) -P 2m-2 (z) coincides with the solution h(z), analytic at z 1 , of the differential equation

m-1 p=0 α p z m-1-p (z 2 + 2az + 1) p h (p) (z) = H(z).
Note that the function h(z) can be explicitely computed as a series expansion around z 1 , see Lemma 4.1. We know that h 1 (z) vanishes at the origin with order m -1, hence it can be written as

h 1 (z) = z m-1 ( h(z) + P m-1 (z)),
where h(z) is known and P m-1 (z) is some polynomial of degree at most m -1. From the link between f (z) and h 1 (z) follows that

z m-1 f (z) = z m-1 ( h(z) + h(1/z)) + Q 2m-2 (z), where Q 2m-2 (z) = z m-1 ( P m-1 (z) + P m-1 (1/z))
is a polynomial of degree at most 2m -2. From the global relation (4.4) and the fact that h(z) is known, we derive that the following integral is also known for µ ∈ A,

S(µ) 2iπ ˆT Q 2m-2 (z) (1 -ϕ(µ)z) m (z -µ) m dz = S(µ) (m -1)! Q 2m-2 (z) (1 -ϕ(µ)z) m (m-1) (µ) = 1 (m -1)! m-1 p=0 α p µ m-1-p (µ 2 + 2aµ + 1) p Q (p) 2m-2 (µ)
Consequently, Q 2m-2 (z) solves a differential equation of the type (4.9) with a second member which is necessarily a polynomial, hence analytic at z 1 . From Lemma 4.1, we know that Q 2m-2 is uniquely determined and can be recovered from this equation. Therefore, P m-1 (z), h 1 (z), g 1 (z), and finally f (z) are also uniquely recovered.

α p z m-1-p (z 2 + 2az + 1) p h (p) (z) = H(z), (4.9)
which is analytic in a neighborhood of z 1 . Moreover, if H(z) = P 2m-2 (z) is a polynomial of degree at most 2m -2, then the solution h(z) exists and is also a polynomial of degree at most 2m -2.

Proof. The two roots z 1 and z 2 of z 2 + 2az + 1 are regular singular points of the differential equation (4.9). We rewrite it in a neighborhood of

z 1 = -a + √ a 2 -1 ∈ D, leading to m-1 p=0 α p (z + z 1 ) m-1-p (z + β) p z p h (p) (z) = H(z + z 1 ), where β = 2 √ a 2 -1 and h(z) = h(z + z 1 ). Denoting by ∞ k=0
a k z k the series expansion of h in a neighborhood of 0, we have

z p h (p) (z) = ∞ k=p k . . . (k -p + 1)a k z k and therefore m-1 p=0 α p (z + z 1 ) m-1-p (z + β) p ∞ k=p k . . . (k -p + 1)a k z k = H(z + z 1 ). (4.10)
We compute the coefficient of z n in the left-hand side of (4.10) which equals

1 n! m-1 p=0 α p (z + z 1 ) m-1-p (z + β) p ∞ k=p a k k! (k -p)! z k (n) (0)
.

By Leibniz rule we get

1 n! m-1 p=0 α p k 1 +k 2 +k 3 =n 1 k 1 !k 2 !k 3 ! (m -1 -p)! (m -1 -p -k 1 )! z m-1-p-k 1 1 p! (p -k 2 )! β p-k 2 a k 3 (k 3 !) 2 (k 3 -p)! ,
where we adopt the convention that 1/q! = 0 if q is a negative integer. Taking into account the definition of α p , this simplifies to Interchanging the two last sums, we obtain

m-1 =0 a n-β n , n ≥ 0,
where a n = 0, n < 0, and we have set, for ∈ {0, . . . , m -1}, (-β/z 1 ) k (2m -k -2)! (m -1 -k)!k!(n -k)! = 0, since β/z 1 < 0. The nonvanishing of the coefficients β n 0 implies that each coefficient a n , n > 0, is uniquely determined whence the first assertion of the lemma.

β n := (n -)! n!β m-1 k=0 (β/z 1 ) k (m -1 -k)!(k -)!
Next, in order to prove the vanishing of all coefficients a n , n ≥ 2m -1, when the right-hand side is a polynomial of degree at most 2m -2, it is sufficient to show that z m+1 dz = 2iπ m! ((z -µ) m (1 -ϕ(µ)z) m g 1 (z)) (m) (0), (4.14) so that, as before, (4.13) can be rewritten as Φ 1 (µ) -Φ 1 (ϕ(µ)). Since this expression only involves the derivatives of g 1 at 0 of order up to m, it does not contain enough information to reconstruct g 1 (z) and a fortiori the function f (z).

β n N +1 = • • • = β n m-1 =

) see ( 3 .

 3 [START_REF] Henrici | On the domain of regularity of generalized axially symmetric potentials[END_REF])-(3.23) for a precise definition of J(z, k). The path of integration in the above integral is the subarc from z r to k on C a . It lies in {Im z ≥ Im z r } when Im k ≥ Im z r and in {Im z ≤ Im z r } when Im k ≤ Im z r . The definitions of J 0 and the differential form W involve the square root λ

Figure 1 :

 1 Figure 1: The paths of integration γ 1 and γ 2 from z r to z, respectively along C up and C low for the definition (3.2) of φ(z, k).

Figure 2 :Figure 3 :

 23 Figure 2: Orientation of contours and path of integration γ 1 (solid line) for the definition (3.14) of φ(z, k) on the upper sheet S z,1 . The path starts at z r , goes along the upper part C up of the circle C a to z l and then follows the horizontal segment up to z. The choice of the determination of the square root λ at the initial point z r of the path γ 1 depends on the region which contains k.

  1) allows one to recover the Neumann data u n from the Dirichlet data u t . 4.1 Case of an even negative integer coefficient α Proof of Theorem 1.5 when α = -2(m -1), m ∈ N * . Equivalently to recovering the Neumann data u n , we may recover the function f (z) := (x + a)u n (z + a), z ∈ T, (4.2)

Lemma 4 . 1 .

 41 Let H(z) be a function analytic in a neighborhood of z 1 . There exists at most one solution h(z) of the differential equation m-1 p=0

(- 1 ) 1 p k 1 3 . 1 p k 1 +k 2 =z -k 1 1 β

 1113121 p (2m -p -2)! β z +k 2 +k 3 =n k 3 !z -k 1 1 β -k 2 k 1 !k 2 !(m -1 -p -k 1 )!(p -k 2 )!(k 3 -p)! a kLet us write k 3 = n -where is a new index. Then k 1 + k 2 = in the previous sum.Since k 1 ≤ m -p -1 and k 2 ≤ p, we deduce that takes its values in the set {0, . . . , m--k 2 k 1 !k 2 !(m -1 -p -k 1 )!(p -k 2 )! , or writing k 1 = k -p and k 2 = p + -k in terms of a new index k, z 1 ) k (k -p)!( + p -k)!(m -1 -k)!(k -)! .

1 p=0(- 1 )

 11 p (2m -p -2)! (k -p)!(n --p)!( + p -k)! . (4.11) If = 0, the second sum in (4.11) equals mp (2m -p -2)! (k -p)!(n -p)!(p -k)! = (-1) k (2m -k -

1 p=0(- 1 )

 11 0, N ∈ {0 . . . , m -2}, n = 2m -1 + N. (4.12)Let N be fixed in {0 . . . , m -2} and ∈ {N + 1, . . . , m -1}. From our convention on the factorials of negative integers, the second sum in (4.11) actually simplifies tok p=k-(-1) p (2m -p -2)! (k -p)!(n --p)!( + p -k)! = (-1) k- p=0 (-1) p (2m -k + -p -2)! ( -p)!p!(n -k -p)! = (-1) k- p=0 (2m -k + -p -2)! (2m -1 + N -k -p)! (-1) p p!( -p)! .The last expression equals the derivative of order -(N + 1) of the polynomial (-1) N +1- p=0x 2m-k+ -p-2 p!( -p)! taken at x = -1. Thanks to the binomial identity, this polynomial rewrites as(-1) N +1-x 2m-k-2 (1 + x) / !.We then get mp (2m -p -2)! (k -p)!(n --p)!( + p -k)! = 0 which shows (4.12).

4. 2

 2 Case of an even positive integer coefficient αHere we assume α = 2(m + 1), m ∈ N. First, let us check that the global relation (4.1) does not allow the reconstruction of u n in this case. Indeed, integrating on T instead of C a , the part involving u n in (4.1) now equalsˆT (z -(k -a)) m (z + 1/(k + a)) m f (z) z m+1 dz (4.13)where f is defined by (4.2) and k ∈ C \ (D a ∪ D -a ). With the definitions of g 1 , µ, ϕ and A as in the previous case, we let, for µ ∈ A, Φ 1 (µ) := ˆT (z -µ) m (1 -ϕ(µ)z) m g 1 (z)

  Proof of Theorem 1.5 when α = 2(m + 1), m ∈ N. Expressing the closed differential form (2.6) in terms of the directional derivatives, thanks to (3.7), and integrating along C a , we now obtain ˆCa((y + ik)xu t (z) -x 2 u n (z) -(2m + 1)(zz + iyk)u(z))x 2m (k -z) m+1 (k + z) m+1 ds = 0, for k ∈ C \ (D a ∪ D -a). Since we assume that u and u t are known on the boundary C a of the domain, our problem is to recover u n from the knowledge, for k∈ C \ (D a ∪ D -a ), of the integral ˆCa x 2m+2 u n (z) (k -z) m+1 (k + z) m+1 ds,which is precisely what was done in the previous case (except that the definition (4.2) of the function f there has to be replaced with f (z) = (x + a) 2m+2 u n (z + a)).

  Assume α = -2(m -1), m ∈ N * , and u t is a given function in L 2 (C a ). Let u n be a function in L 2 (C a ) such that the global relation (1.14) holds true. Then, u n is unique and can be explicitly recovered from that relation. A similar statement holds when α = 2(m + 1), m ∈ N. In that case, (1.14) does not allow the reconstruction of u n . A relation that works is obtained by integrating another differential form, see(2.6). Making use of the equation conjugate to(1.6), that is changing α into -α, we derive the converse reconstruction of the Dirichlet data u t from the Neumann data u n , still when α ∈ 2Z is an even coefficient. In Theorem 1.5, it is actually not necessary to assume smoothness of the functions u t and u n on the boundary circle C a . Hence, we only assume these functions to be L 2 on C a .

	14)
	It is sometimes conjectured that there always exists at least one global relation (i.e. a
	Lax pair) that allows for the recovering of one type of boundary values from the other
	one. In the case of an even coefficient α, we show that, indeed, this correspondance can
	be performed explicitly from such a relation.
	Theorem 1.5. Remark 1.6.

  Proposition 2.1. The function u satisfies (2.1) in a bounded, simply connected domain Ω of H if and only if the differential form


Céline Teulière
email: celine.teuliere@univ-bpclermont.fr

Eric Marchand
email: marchand@irisa.fr

A dense and direct approach to visual servoing using depth maps

Keywords: dense sensor-based control, depth map, visual servoing

This paper presents a novel 3D servoing approach using dense depth maps to perform robotic tasks. With respect to positionbased approaches, our method does not require the estimation of the 3D pose (direct), nor the extraction and matching of 3D features (dense) and only requires dense depth maps provided by 3D sensors. Our approach has been validated in various servoing experiments using the depth information from a low cost RGB-D sensor. Positioning tasks are properly achieved despite noisy measurements, even when partial occlusions or scene modifications occur. We also show that, in cases where a reference depth map cannot be easily available, synthetic ones generated with a rendering engine still lead to satisfactory positioning performances. Application of the approach to the navigation of a mobile robot is also demonstrated.

I. INTRODUCTION

Most of the robotic positioning tasks are still achieved today by estimating first the relative pose between the robot and the scene or the object of interest, and then using a pose-based control scheme, as initially proposed in [START_REF] Wilson | Relative End-Effector Control Using Cartesian Position Based Visual Servoing[END_REF]. However, the pose estimation problem itself is complex in its general formulation. Also known as the 3D localization problem [START_REF] Lepetit | Monocular Model-Based 3D Tracking of Rigid Objects: A Survey[END_REF], this problem has been widely investigated by the computer vision community but remains non-trivial for vision sensors alone, in particular in low-textured environments. Using range data, a range flow formulation has been proposed [START_REF] Horn | Rigid body motion from range image sequences[END_REF] [START_REF] Gharavi | 3-D Motion Estimation Using Range Data[END_REF] to estimate the 3D pose of a mobile robot. Alternatively, the alignment of successive 3D point clouds using ICP [START_REF] Besl | A method for registration of 3-D shapes[END_REF], [START_REF] Chen | Object modeling by registration of multiple range images[END_REF] has become a very popular method. Many variants have been proposed in the literature [START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF] and the development of the so-called RGB-D cameras attracted a lot of attention on these methods in the recent years [START_REF] Swadzba | A comprehensive system for 3D modeling from range images acquired from a 3D ToF sensor[END_REF], [START_REF] May | Three-dimensional mapping with time-of-flight cameras[END_REF], [START_REF] Henry | RGB-D Mapping: Using depth cameras for dense 3D modeling of indoor environments[END_REF], [START_REF] Ra | KinectFusion: Real-time dense surface mapping and tracking[END_REF]. Some work also consider crude global alignment using a global voting scheme within a transformed space [START_REF] Makadia | Fully automatic registration of 3d point clouds[END_REF].

In this paper, we propose to perform robotic tasks without reconstructing the full 3D pose between the robot and its environment, but using a sensor-based servoing scheme, the considered data being directly the depth map obtained from a range sensor. Our approach is thus related to other sensor-based methods, such as image-based visual servoing (IBVS) [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF], where a robotic task is expressed directly as the regulation of a visual error. In IBVS, the visual error is usually defined as the difference between a current and a desired set of geometric features (points, straigth lines, etc.) selected from the image, to control the desired degrees of freedom. Therefore, IBVS schemes usually require the extraction of visual features from image measurements, and their matching in successive frames. However, those steps, based on image processing techniques, are often considered as the bottleneck of visual servoing methods. In the tracking literature, dense approaches that do not require matching have already been proposed, based on the "brightness constancy constraint" [START_REF] Hager | Efficient region tracking with parametric models of geometry and illumination[END_REF], [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF], [START_REF] Silveira | Real-time visual tracking under arbitrary illumination changes[END_REF]. In [START_REF] Corso | Direct plane tracking in stereo images for mobile navigation[END_REF] this constraint is used in a stereo system to track planes and is applied to mobile robot localization. Recently, some visual servoing work also proposed to use all the image directly, without any extraction or matching step, by minimizing the difference between the current image and a reference image. This approach is referenced as photometric visual servoing [START_REF] Collewet | Visual servoing set free from image processing[END_REF] [START_REF] Collewet | Photometric visual servoing[END_REF]. However, luminance-based approaches are not always applicable since they require texture, and stable lighting condition, or an accurate knowledge of the materials and light sources to model image formation [START_REF] Silveira | Unified direct visual tracking of rigid and deformable surfaces under generic illumination changes in grayscale and color images[END_REF].

In our work we propose to use the dense depth map obtained from a range sensor as a visual feature for positioning task wrt. non planar scenes, without any feature extraction or matching step. This is a major difference with respect to approaches such as [START_REF] Martinet | Vision based control law using 3D visual features[END_REF] where 3D points have to be matched. We derive a new control law for a robot positioning or navigation using this feature directly. Our approach is thus both direct (without any 3D pose estimation) and dense (without feature extraction or matching). A first version of this work has been presented in [START_REF] Teulière | Direct 3D servoing using dense depth maps[END_REF]. We provide here an extended version, with new experimental results underlying the strong potential of this approach. To the best of our knowledge, this is the first work proposing such a dense depth-based visual servoing.

II. DIRECT DENSE DEPTH MAP SERVOING

This section presents the heart of our approach, i.e. how to control a robot using dense depth maps. We first introduce what we call a depth map and what it means to use it as feature (Section II-A). Then we derive the fundamental equations necessary to compute our control law (Sections II-B and II-C).

A. Depth map sensing

There are multiple technologies of sensors capable of providing depth (or range) information. Most range sensors without contact are active, and based on the time of flight (ToF) principle: the idea is to send waves of known velocity and measure the time it takes them to go from the sensor and come back after reflection on the scene. This can be achieved by sending light pulses. Another approach consists in using a modulated signal and measuring the phase shift. In each case, the depth information is derived knowing the velocity of the sent signal (eg: Laser scans, sonars, radars, ToF or RGB-D cameras). Another existing technology for active range sensing is based on structured light: known patterns (stripes, dots, ...) are projected onto the scene and the depth information is deduced from their deformation. This technology is used for instance in the recent Microsoft Kinect or Asus Xtion pro devices, based on PrimeSense technology [START_REF] Freedman | Depth mapping using projected patterns[END_REF]. Depth can also be measured with passive sensors such as cameras: by matching image features in two different views of a calibrated stereo rig, depth can be computed from geometry. The depth information is sparse when a finite set of features are matched, but dense depth maps can also be obtained [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF].

In the following, we consider a range sensor capable of providing dense depth maps. Without loss of generality, the range measurements are expressed in sensor centered cartesian coordinates. We also consider that the depth map is represented according to a perspective projection model (see Figure 1). This is a natural choice for any range sensor based on perspective cameras (stereo pairs, or structured light such as the Kinect sensor used in our experiments). It is also very general since any dense depth map coming from other sensors (laser, radar, etc.) can be converted with such a perspective projection. Formally, we denote by Z(x, y, t) the depth at time t of the 3D point of coordinates (X, Y, Z) in the sensor frame, with X = xZ and Y = yZ, (x, y) being the metric image coordinates. Figure 2 gives an example of depth map obtained from Microsoft Kinect RGB-D sensor, where the depth values have been scaled to greyscale levels. White pixels correspond to unavailable depth values, i.e. pixels where the sensor could not compute any depth information. Note also that for better visualisation purpose, we applied histogram equalization The next section shows how such dense depth maps can be used to control a robot.

B. Modeling

Let us consider that a robot end effector is equipped with a range sensor (Figure 3). We express a positioning task as the regulation of the feature Z to a desired value Z * . Here, Z = (Z1, ..., ZN ) is a vector containing the N depth values corresponding to the current dense depth map. The desired value Z * thus corresponds to a reference depth map acquired at the desired robot pose.

Therefore, the control law to design aims at regulating the following error to zero:

e = Z -Z * =     . . . Zi -Z * i . . .     (1)
An illustration of such an error is given in Figure 4.

In analogy with the visual servoing framework [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF] we denote by L Z the interaction matrix associated to the feature Z, and characterized by the relation:

∂Z ∂t = L Z v (2)
where ∂Z ∂t is the temporal variation of the depth and v = (v, ω) is the sensor instantaneous velocity, with v the translational velocity, and ω the rotational velocity.

We now derive the expression of the matrix L Z which will be required in the control law (Section II-C). In the following, we consider the continuous formulation of the depth map as a surface Z(x, y, t). Assuming that the scene is rigid and the surface Z(x, y, t) is smooth1 , taking its full derivative leads to:

Ż = dZ dt = ∂Z ∂x ẋ + ∂Z ∂y ẏ + ∂Z ∂t , (3)
where ( ẋ, ẏ) is the 2D velocity of the image point (x, y). Equation ( 3) is known as the range flow constraint equation [START_REF] Yamamoto | Direct estimation of range flow on deformable shape from a video rate range camera[END_REF] or elevation rate constraint equation [START_REF] Horn | Rigid body motion from range image sequences[END_REF]. It is very similar to the brightness change constraint equation that is used in the computation of optical flow [START_REF] Horn | Determining optical flow[END_REF] and used in direct photometric visual servoing methods [START_REF] Collewet | Photometric visual servoing[END_REF]. The main difference is that in the luminance case, an additional assumption is made to constrain the brightness time derivative to be zero.

From equation ( 3), the temporal variation of the depth is immediately deduced:

∂Z ∂t = Ż -A ẋ -B ẏ, (4)
where A = ∂Z ∂x and B = ∂Z ∂y . Therefore, the interaction matrix LZ related to one depth value is expressed by:

LZ = LP Z -ALx -BLy. (5)
The matrices Lx, Ly defined such that ẋ = Lxv and ẏ = Lyv are the well-known interaction matrices of image point coordinates, given by [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]:

Lx = -1 Z 0 x Z xy -(1 + x 2 ) y (6) Ly = 0 -1 Z y Z 1 + y 2 -xy -x , (7)
and LP Z is the interaction matrix related to the coordinate Z of a 3D point, such that Ż = LP Z v. It is given by [3][20]:

LP Z = 0 0 -1 -yZ xZ 0 . (8)
Finally, replacing ( 6), ( 7) and ( 8) in ( 5), we get:

LZ = A Z B Z -Z+xA+yB Z Zw x Zw y Zw z , (9)
where Zw x = -yZ-xyA-(1+y 2 )B, Zw y = xZ+(1+x 2 )A+xyB and Zw z = xB -yA. Note that this expression underlines one of the main differences between our approach and sparse 3D approaches [START_REF] Martinet | Vision based control law using 3D visual features[END_REF], as discussed in Appendix A. The full interaction matrix L Z of size N × 6 corresponding to the entire depth map is thus the stack of the 1 × 6 matrices LZ i :

L Z =    LZ 1 . . . LZ N    . (10)

C. Control law

We consider the following control law:

v = -λL + Z (Z -Z * ) ( 11
)
where λ is a scalar gain parameter and L + Z denotes the pseudo-inverse of L Z defined by

L + Z = (L Z L Z ) -1 L Z .
Note that exactly the same demonstration as for IBVS as given in [START_REF] Chaumette | Visual servoing and visual tracking[END_REF] demonstrating the local asymptotic stability (LAS) of this control law under the condition that L + Z L Z > 0 in the neighborhood of Z = Z * . This condition is ensured as soon as L Z is of full rank 6 since L + Z L Z = I6 in that case. The interaction matrix given by [START_REF] Freedman | Depth mapping using projected patterns[END_REF] has to be of full rank for the system to be stable. For instance, a fully planar scene would lead to a rank 3 matrix where only 3 degrees of freedom (dof) could be controlled to form a plane-to-plane virtual linkage. The interpretation of this is that an infinite number of poses would lead to the same depth map perception. However the full-rank assumption is easily fulfilled in real world scenarios where depth variations are discriminative enough to avoid indetermination. Typically, observing points from 3 non parallel planes is sufficient.

III. PRACTICAL ISSUES AND ROBUSTNESS IMPROVEMENTS

In the previous section, we presented our depth map based servoing method. When testing it, we found that this method was efficient in simulation sequences, with perfect data, but we had to face some practical issues in real conditions, in particular, in our case, using a Kinect sensor. This section presents the modifications we had to undertake in order to improve the robustness of the servoing task with respect to noisy and incomplete measurements (Section III-A) and to scene perturbations and occlusions (Section III-B).

A. Noisy and incomplete measurements

As illustrated in Figure 2-b the depth map acquired by a Kinect sensor is noisy and incomplete. In practice, we only considered the pixels for which a depth value was available both in the reference Z * and the current Z depth maps. This means that the number N of depth values in Z and [START_REF] Freedman | Depth mapping using projected patterns[END_REF], is inferior to the size of the depth map.

In the experiments presented in this paper, about 80% of the total number of pixels could typically be used.

In addition, we reduced the noise by applying a standard 3 × 3 Gaussian filter on the depth maps, the convolution being computed only with the valid neighbors.

Similarly, the spatial gradient was computed using a standard 3×3 derivative kernel taking into account the valid neighbors only.

B. Occlusions and scene modifications

Another issue to take into account is the possibility of partial occlusions or scene modifications during the servoing process. To reduce the effect of such events on the task achievement, we use robust M-estimation [START_REF] Huber | Robust Statistics[END_REF]. We thus introduce a modification of our task objective (1) allowing uncertain measures to be less likely considered or in some cases completely rejected. The new task error is given by [START_REF] Comport | Statistically robust 2-D visual servoing[END_REF]:

e = D(Z -Z * ) ( 12
)
where D is a diagonal weighting matrix: D = diag (w1, ..., wN ). The new control law thus becomes:

v = -λ(DL Z ) + D(Z -Z * ), (13)
The weights wi in D are computed using Tukey's estimator. The reader can refer to [START_REF] Huber | Robust Statistics[END_REF] for details on M-estimation.

IV. EXPERIMENTAL RESULTS

In this section we first provide the experimental validation of our approach for positioning tasks (Section IV-A). Then its application to the navigation of a mobile robot is proposed (Section IV-B).

A. Positioning tasks

In our positioning tasks, a Kinect sensor has been mounted on a ADEPT Viper robot (see Figure 3). In each experiment, the task is expressed as the minimization of the error (12) between a fixed desired depth map and the current one. The control law ( 13) is computed with a fixed gain λ = 2.5. The depth maps are acquired using the LibFreenect2 driver through the ViSP library [START_REF] Marchand | ViSP for visual servoing: A generic software platform with a wide class of robot control skills[END_REF], with a resolution of 320 × 240 pixels.

In terms of computation, each iteration requires the computation of 3 × 3 gradients in each (non-discarded) pixel which is very fast. The most costly step is to fill the N × 6 interaction matrix L Z to compute the control law. Note that L Z L Z used to compute the pseudo-inverse is a 6 × 6 matrix which is very fast to invert. Without any specific optimization, the code runs in about 60 ms per frame on a standard laptop. The method is suitable for real-time experiments.

1) Using a synthetic depth image to define the desired position: In the first experiment, we consider the case where the desired position is defined in a simulation environment using a model of the scene. This kind of approach can typically be beneficial in applications where a 3D model is known but one wants to define different tasks without the need for depth maps acquired in situ. In that case we render the desired depth image from a 3D model of the scene instead of using one acquired from the sensor.

For this experiment, we built a 3D model of the scene using the ReconstructMe3 software. We then defined the 3D pose we wanted the robot to reach, with respect to this model. The desired depth map corresponding to this pose was rendered using Ogre3D4 using the actual depth camera calibration parameters. In this case the depth sensor thus needs to be calibrated. The positioning task consists for the robot in minimizing the error between this rendered depth map and the current one acquired by the Kinect sensor.

The results of this experiment are shown in the Figures 6 and7. The first image (1-a) shows the RGB view provided by the Kinect for the initial position. This image is never used in the control scheme and is only given here for a better understanding of the setup. The depth maps are shown in the second row, and the last row gives the corresponding error, i.e. the difference between the desired and the current depth maps, unavailable data being discarded as explained in III-A. The difference images are scaled so that a plain grey frame (3-b) corresponds to a null error, and thus to the good achievement of the task. Figure 6 (3-a) gives a vizualization of the error in the initial position.

The corresponding quantitative values for the task error, the 3D positioning errors and the velocities are given in Figure 7. Figure 7 (b) shows the repartition of the depth errors in the initial frame of figure 6 (3-a). The maximum error is about 50 cm in depth. The peak at 0 corresponds mostly to white pixels of Figure 6 (1-b) that are not considered.

Note that one iteration corresponds to one execution of the control loop, that is the computation of (13) for the current image and the transmission of this velocity command to the robot. Figures (c) and (d) show that with an initial error of 15cm in translation and 10 • in rotation, the positioning task is properly performed, as indicated by the low residuals. Here the final accuracy depends on the quality of the model. To evaluate the final accuracy of the system itself we thus consider in the next experiment a desired depth map acquired at the desired position.

2) Robustness to occlusions and scene changes: In the second experiment, we evaluate the robustness of our approach with respect to partial occlusions or modifications of the observed scene. First, the desired depth map is acquired at the desired position, then the robot is moved to the initial state in which the servoing is launched. The goal here is to regulate the error between this reference depth map and the current depth map. The initial scene is illustrated in Figure 9 (1-a). During the task achievement, someone entered the sensor field, removed an object and put it back several times. Some selected frames of this sequence are shown in Figure 9. The full video of this sequence is provided as supplementary material. The initial and final positions are illustrated in the first and last columns, while columns (b) and (c) show examples of occlusions. Note that at the end of the sequence the white bear has been completely removed from the scene, and the final depth map (Figure 9 (2-d)) is thus different from the desired one (Figure 8 (b)). This difference appears in the final difference image (Figure 9 (4-d)) and the task error function (Figure 8 (a)). However, despite the scene modifications and occlusions, the positioning task is successfully achieved, as shown by the convergence of the positioning errors in Figure 8 (c). The robustness of our control scheme to perturbations is the result of the use of M-estimation (see III-B). The effect of M-estimation is illustrated on the third row of Figure 9, where the relative weights of each data in equation ( 12) are represented. Black pixels correspond to rejected values and brightest ones to inliers. Figures 8 (c) show that with an initial error of 17cm in translation and 20 • in rotation, the positioning task is properly achieved with a remaining error of less than 3mm in translation and 0.4 • in rotation. Given the low depth resolution of the sensor and its noisy measurements, this corresponds to a good achievement of the task.

Finally, note that in this scene the smoothness assumption was not verified everywhere since large depth discontinuities exist at the border of the objects, for example between the table and the floor. This experiment thus shows that the method is successful beyond its initial assumption. This is due to the fact that the points corresponding to discontinuities are a minority and are detected as outliers by the M-estimator. 3) Experimental analysis of the convergence domain: We want to underline that the above experiments are examples from a large number of experiments that we performed using different initial poses and scenes. As for IBVS approaches, only the local asymptotic stability can be demonstrated (see [START_REF] Chaumette | Visual servoing and visual tracking[END_REF]). Thus we can only assess the performances in terms of convergence domain from experiments. In order to empirically evaluate the convergence domain we run multiple positioning tests from different initial positions and recorded the convergence succes and failure. One test is deemed to have converged if both the norm of the transation error vector and the norm of the rotation error vector get smaller than 1 cm and 1 • respectively, in less than 500 iterations. Those tests were performed using the 3D model of Figure 5 in simulation, without adding occlusions. Simulation allows us to handle exhaustive testing with hundreds of different initial positions. The initial poses are chosen so that the camera center is placed on a regular 3D grid centered on the desired pose in x and y, and with z varying from 0.3 m to 4.7 m where the desired z camera position corresponds to 1.5 m. The orientation is set so that the desired and initial depth maps overlap. This setting leads to large variations of x-axis and y-axis rotations, from -60 • to 60 • . We also considered rotation around z-axis, by running one full set of simulations with 0 • z-axis rotation and another one with 30 • zaxis rotation. Figure 10 shows the resulting convergence domains, which as can be seen are considerable. Note that joint limits are not considered in this simulation test. From Figure 10 we can see that the convergence domain has an approximate radius of more than 1 meter along the x and y directions in this setting, and even more in the z axis where convergence can be obtained from more than 2 m above the desired position. The methods handles indeed easily a large initial error on the depth axis which generates a large velocity component on this axis to compensate for the important depth map error. Note that the convergence domain is scene-dependant, and large structural elements with smoothly varying depth (planes or large rocks in this example) will generally lead to a larger convergence domain than scenes with high frequency depth variations.

B. Non honolomic robot navigation using depth map memory

In the experiments above, we considered positioning tasks on a 6 dof robot, in various conditions. In this part, we propose to apply our depth-based approach to a navigation task on a wheeled robot (Pioneer P3-DX). We assume here that we have stored a sequence of depth maps acquired during a manual navigation stage. These depth maps can be considered as a sensory memory that is then used for the autonomous navigation. The navigation task is thus defined as a succession of positioning tasks using the successive stored data as desired depth maps (see Figure 11). In this case, each positioning subtask is performed in a similar manner as in the previous experiments, but controlling 2 dof only: the forward translational motion and the in-plane rotation. The switch from one reference depth map to the next is based on a simple threshold on the error decrease.

Figure 12 shows some samples from our navigation sequence. The first row represents the robot during the learning phase, that is when it was manually controlled and acquiring the reference depth maps. The second row gives the autonomous navigation results. Although no ground truth measurement was available for this task, this figure shows that the robot closely follows the path of the learning phase, using the stored maps as references.

Figure 13 provides typical depth maps from this navigation sequence, along with the corresponding errors when the reference frame is changed. Note that since the robot is non-holonomic, one cannot ensure its convergence to the 3D position corresponding to the desired depth map, which explains that the error images are not as good as for the previous 6 dof positioning tasks. Note also that one possible indetermination can occurr in the case of a long corridor with no Sensory path Fig. 11. A sequence of depth maps has been acquired in a manual navigation step. The navigation task is then defined using each depth map from the memory as an intermediate desired depth map for the controller.

door to mark a depth variation, since several different poses would lead to the same depth map (the interaction matrix would then be rank deficient). This however would not be an issue for tasks such as moving forward while being centered in a corridor. This experiment shows the feasibility of the approach for indoor navigation.

V. CONCLUSIONS

We have demonstrated that it is possible to use a dense depth map directly to control robot motion. The goal position can be defined by a single depth map either directly acquired from that position or synthetically rendered. The main advantage of our approach is that it does not require any pose estimation, feature extraction or matching step. Moreover, when the depth map is obtained from an active sensor, the resulting approach is not sensitive to illumination changes as photometric approaches can be. Some limitations can appear with the use of active sensors such as Kinect RGB-D camera, in particular the noise and the absence of some measurements. We show however that those issues can be overcome thanks to the use of M-estimators and basic image pre-processing. APPENDIX Depth information has already been used in position-based visual servoing. For example, [START_REF] Martinet | Vision based control law using 3D visual features[END_REF] proposed to use the 3D coordinates (X, Y, Z) of a set of 3D points as features to be regulated in a proportional control law. In other words, the positioning task was expressed as the regulation of the feature P = (X1, Y1, Z1, ..., XN , YN , ZN ) to a reference feature P * = (X * 1 , Y * 1 , Z * 1 , ..., X * N , Y * N , Z * N ) corresponding to the 3D coordinates of the set of points at the desired robot position. The interaction matrix related to a single 3D point is then given by [3][20]:

LP =   -1 0 0 0 -Z Y 0 -1 0 Z 0 -X 0 0 -1 -Y X 0   . ( 14
)
At first sight, the depth components of this kind of 3D feature (X1, Y1, Z1, ..., XN , YN , ZN ) could seem very close to the vector formulation Z = (Z1, ..., ZN ) that we defined in Section II-C. However, a key difference with respect to our approach is that [START_REF] Martinet | Vision based control law using 3D visual features[END_REF] uses a sparse set of 3D features. Consequently, in [START_REF] Martinet | Vision based control law using 3D visual features[END_REF] a matching step is required to determine the feature values through the sequence, and the range flow equation (3), based on a smoothness assumption, does not hold in the sparse case. On the contrary, one of the key advantages of the method we propose, is that it does not require any feature extraction nor matching step and uses directly the dense depth information from the range sensor thanks to the range flow equation. That is why the interaction matrix related to the depth map is given by ( 9) while the interaction matrix related to the depth of point is given by the last row of ( 14).

Fig. 1 .

 1 Fig. 1. Sensor frame representation.

Fig. 2 .

 2 Fig. 2. Example of static scene (a) and corresponding depth map representation (b) acquired from Microsoft Kinect sensor. The darkest pixels correspond to the smallest depths. White pixels correspond to unavailable data.

Fig. 3 .

 3 Fig. 3. ADEPT Viper robotic system equipped with a Microsoft Kinect sensor.

Fig. 4 .

 4 Fig. 4. The task error is the difference of depth maps Z -Z * .

Figure 6 (

 6 Figure6(1-b) shows an example of such a rendered image for the model of Figure5. The white points in this image correspond to unmodeled areas, for which the rendering gives an infinite depth value. They are treated as unavailable data, as for Kinect depth maps, and are excluded from the feature set (see Section III).

Fig. 5 .

 5 Fig. 5. 3D model used for generating the desired depth map in Figure 6.

Fig. 6 .

 6 Fig. 6. experiment. First column corresponds to the initial position. The RGB view from the Kinect (1-a) is not used in the algorithm. (2-a) Initial depth map, where white parts correspond to unavailable data. (3-a) Difference between the initial and desired depth maps. Second column corresponds to the end of the motion. The desired depth map (1-b) was rendered using a 3D model of the scene. (2-b) shows the final depth map, which minimizes the difference (3-b) with the desired one (1-b).

Fig. 7 .

 7 Fig. 7. First experiment. (a) Task error, (b) histogram of depth errors (m) corresponding to the initial frame, (c) translational part of positioning error, (d) rotational part of positioning error, (e) translational velocities, (f) rotational velocities.

Figure 9 (

 9 3-b), (3-c), and (3-d) show that the perturbations are correctly detected since the corresponding pixels are given a smaller weight.

Fig. 8 .

 8 Fig. 8. Second experiment. (a) Task error, (b) desired depth map, (c) positioning errors, (d) velocities.

Fig. 10 .

 10 Fig. 10. Experimental evaluation of the convergence domain. Result are given for zero z-axis rotation in column (1) and 30 • z-axis rotation in column (2).The first row shows for both cases the convex hull of the points from where the system successfully converged. Rows (b) and (c) shows slices at the desired y and z planes respectively. The poses that converged are shown by green circles and those that diverged by red crosses.

Fig. 9 .

 9 Fig. 9. Selected frames of the second experiment. Columns correspond to frames 1, 15, 69 and final frame respectively. Those frames illustrate occlusions and object removal (1-b) (1-c) (1-d). The first row gives the RGB view from the Kinect, which is not used in the algorithm but shows the setup. The depth maps are represented in the second row. The images of the third row represent the weights of each pixel in the M-estimation. Black pixels are discarded. Frames (3-b) (3-c) (3-d) show that occluded areas are given a very low weight. Fourth row: difference between the initial and desired depth maps.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Extracted frames from the navigation experiment. The first row corresponds to the manual navigation step where the depth maps are memorised. The second row shows the autonomous navigation using this sensory memory.

The points in the surface where this assumption do not hold will be discarded by the M-estimation process presented in section III.

http://openkinect.org/

http://reconstructme.net

http://www.ogre3d.org/

ACKNOWLEDGMENT

This work was supported by ReV-TV FUI. The authors would like to thank Franc ¸ois Chaumette for helpful discussions.