Roman Novikov
email: novikov@math.univ-nantes.fr

Vladimir Sharafutdinov

On problem of polarization tomography, I

The polarization tomography problem consists of recovering a matrix function f from the fundamental matrix of the equation Dη/dt = π γ f η known for every geodesic γ of a given Riemannian metric. Here π γ is the orthogonal projection onto the hyperplane γ⊥ . The problem arises in optical tomography of slightly anisotropic media. The local uniqueness theorem is proved: a C 1 -small function f can be recovered from the data uniquely up to a natural obstruction. A partial global result is obtained in the case of the Euclidean metric on R 3 .

Introduction

First of all we shortly recall the physical motivation of the problem. See Section 5.1 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF] for a detailed discussion.

We consider propagation of time-harmonic electromagnetic waves of frequency ω in a medium with the zero conductivity, unit magnetic permeability, and the dielectric permeability tensor of the form

ε ij = n 2 δ ij + 1 k χ ij , (1.1)
where k = ω/c is the wave number, c being the light velocity. Here n > 0 is a function of a point x ∈ R 3 , and the tensor χ ij = χ ij (x) determines a small anisotropy of the medium. The smallness is emphasized by the factor 1/k. Equation (1.1) was suggested by Yu. Kravtsov [START_REF] Yu | Quasi-isotropic" approximation of geometric optics[END_REF]. By some physical arguments [START_REF] Yu | Geometric Optics of Inhomogeneous Media[END_REF], the tensor χ must be Hermitian, χ ij = χji .

In the scope of the zero approximation of geometric optics, propagation of electromagnetic waves in such media is described as follows. Exactly as in the background isotropic medium, light rays are geodesics of the Riemannian metric dt 2 = n 2 (x)|dx| 2 ;

(1.2) the electric vector E(x) and magnetic vector H(x) are orthogonal to each other as well as to the ray; and the amplitude A 2 = |E| 2 = |H| 2 satisfies A = C/ √ nJ along a ray, where J is the geometric divergence and the constant C depends on the ray. The only difference between a slightly anisotropic medium and the background isotropic one consists of the wave polarization. The polarization vector η = n -1 A -1 E satisfies the equation (generalized Rytov's law)

Dη dt = i 2n 2 π γ χη (1.3)
along a geodesic ray γ(t). Here t is the arc length of γ in metric (1.2), γ = dγ/dt is the speed vector of γ, π γ is the orthogonal projection onto the plane γ⊥ , and D/dt = γk ∇ k is the covariant derivative along γ in metric (1.2). The right-hand side of (1.3) is understood as follows: π γ and χ are considered as linear operators on T C γ(t) = C 3 , and π γ χη is the result of action of the operator π γ χ on the complex vector η ∈ γ⊥ . Here γ⊥ is the twodimensional complex vector space consisting of complex vectors orthogonal to the real vector γ ∈ R 3 = T γ(t) ⊂ T C γ(t) = C 3 . Introducing the notation

f = i 2n 2 χ, (1.4)
we rewrite (1.3) in the form Dη dt = π γ f η.

(

Observe that f is a skew-Hermitian operator, f ij = -fji .

Let us now consider the inverse problem. Assume a medium under investigation to be contained in a bounded domain D ⊂ R 3 with a smooth boundary. The background isotropic medium is assumed to be known, i.e., metric (1.2) is given. The domain D is assumed to be convex with respect to the metric, i.e., for any two boundary points x 0 , x 1 ∈ ∂D, there exists a unique geodesic γ : [0, 1] → D such that γ(0) = x 0 , γ(1) = x 1 . We consider the inverse problem of determining the anisotropic part χ ij of the dielectric permeability tensor or, equivalently, of determining the tensor f on (1.5). To this end we can fulfill tomographic measurements of the following type. For any unit speed geodesic γ : [0, l] → D between boundary points, we can choose an initial value η 0 = η(0) ∈ γ⊥ (0) of the polarization vector and measure the final value η 1 = η(l) ∈ γ⊥ (l) of the solution to equation (1.5). In other words, we assume the linear operator γ⊥ (0) → γ⊥ (l), η 0 → η 1 to be known for every unit speed geodesic γ : [0, l] → D between boundary points. Instead of (1.5), we will consider the corresponding operator equation

D Ũ(t) dt = f γ(t) Ũ(t), (1.6)
where f γ(t) : γ⊥ (t) → γ⊥ (t) is the restriction of the operator π γ(t) f (γ(t)) to the plane γ⊥ (t), and the solution is considered as a linear operator Ũ(t) : γ⊥ (t) → γ⊥ (t). Equation (1.6) has a unique solution satisfying the initial condition

Ũ (0) = E, (1.7)
where E is the identity operator. Since f γ(t) is a skew-Hermitian operator, the solution Ũ(t) is a unitary operator. The final value of the solution

Φ[f ](γ) = Ũ (l) ∈ GL( γ⊥ (l))
is the data for the inverse problem. Given the function Φ[f ] on the set of unit speed geodesics between boundary points, we have to determine the tensor field f = (f ij (x)) on the domain D.

We consider the inverse problem in a more general setting. Instead of a domain D ⊂ R 3 with metric (1.2), we will consider a compact Riemannian manifold (M, g) of an arbitrary dimension n ≥ 3, and an arbitrary complex tensor field f = (f ij ) on M. In such a setting, equation (1.6) makes sense along a geodesic γ. We will subordinate the manifold (M, g) to some conditions that guarantee smoothness of the data Φ[f ] in the case of a smooth f .

The two-dimensional case of n = 2 is not interesting since f γ and Ũ become scalar functions and the solution to the scalar equation (1.6) is given by an explicit formula in this case. Therefore the inverse problem is reduced to the inversion of the ray transform I on second rank tensor fields, see the remark before Theorem 5.2.1 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF].

Equation (1.6) can be slightly simplified. For a point x ∈ M, let T C

x M be the complexification of the tangent space T x M. Instead of considering the operator Ũ(t) on γ⊥ (t), we define the linear operator

U(t) : T C γ(t) M → T C γ(t) M by U(t)| γ⊥ (t) = Ũ (t), U(t) γ(t) = γ(t).
If Ũ (t) satisfies (1.6)-(1.7), then U(t) solves the initial value problem

DU dt = (π γ f π γ )U, U(0) = E. (1.8)
Equation (1.8) is more handy than (1.6) since all operators participating in (1.8) are defined on the whole of T C γ(t) M. The inverse problem consists of recovering the tensor field f from the data Φ[f ](γ) = U(l) known for every unit speed geodesic γ : [0, l] → M between boundary points. Let us emphasize that the inverse problem is strongly nonlinear, i.e., the data Φ[f ] depends on f in a nonlinear manner.

The three-dimensional case, n = dim M = 3, is of the most importance for applications as we have shown above. On the other hand, the three-dimensional case is mathematically the exceptional one because, for a skew-symmetric f , the solution to the inverse problem is not unique. The non-uniqueness is discussed in Section 4.

The main result of the present article is the local uniqueness theorem: the solution to the inverse problem is unique (up to a natural obstruction in the three-dimensional case) if the tensor field f is C 1 -small. See Theorem 5.1 below for the precise statement.

Our method of investigating the inverse problem is a combination of approaches used in [START_REF] Vertgeim | Integral geometry with matrix weight, and a nonlinear problem of recovering matrices[END_REF] and in Chapter 5 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF]. First of all, following [START_REF] Vertgeim | Integral geometry with matrix weight, and a nonlinear problem of recovering matrices[END_REF], we reduce our nonlinear problem to a linear one as follows. Let f i (i = 1, 2) be two tensor fields and U i (t) be the corresponding solutions to (1.8)

with f = f i . Then u = U -1 1 U 2 -E satisfies Du dt = pπ γ (f 2 -f 1 )π γ q, u(0) = 0, (1.9)
where p = U -1 1 and q = U 2 . We consider p and q as operator-valued weights which are close to the unit operator if f i are C 1 -small. Assuming the weights p and q to be fixed, the solution u(t) to the initial value problem (1.9) depends linearly on f = f 2 -f 1 . We study the linear inverse problem of recovering the tensor field f = f 2 -f 1 from the data F [f ](γ) = u(l) given for all geodesics γ : [0, l] → M between boundary points. In the case of a symmetric tensor field f and of unit weights, this linear problem was considered in Chapter 5 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF]. We will demonstrate that the same approach works in the case of an arbitrary f and of weights close to the unit one.

There is one more opportunity to extract a linear inverse problem from equation (1.8). Indeed, if W (t) = det U(t) is the Wronskian, then the function ϕ(t) = ln W (t) satisfies

dϕ(t) dt = tr (π γ f π γ ).
Therefore, for every unit speed geodesic γ : [0, l] → M between boundary points, the integral

S[f ](γ) = l 0 tr (π γ(t) f (γ(t))π γ(t) )dt (1.10) is expressed through the data Φ[f ] by the formula S[f ](γ) = ln det Φ[f ](γ).
The data S[f ] depends linearly on f . Of course, some information is lost while the data Φ[f ] is replaced with S[f ]. In particular, S[f ] is independent of the skew-symmetric part of f . We finally note that main results of the article are new and nontrivial in the case of M ⊂ R n , n ≥ 3, with the standard Euclidean metric. If a reader is not familiar with the tensor analysis machinery on the tangent bundle of a Riemannian manifold, he/she can first read the article for the latter simplest case.

The article is organized as follows. Section 2 contains some preliminaries concerning Riemannian geometry and tensor analysis. In particular, we define some class of Riemannian manifolds for which the problem can be posed in the most natural way. Instead of considering the ordinary differential equation (1.8) along individual geodesics, we introduce a partial differential equation on the unit tangent bundle and pose an equivalent version of the problem in terms of the latter equation. In Section 3, we consider the corresponding linear problem and prove the uniqueness for weights sufficiently close to the unit in the case of n ≥ 4. Section 4 discusses the three-dimensional case. In Section 5, we check that the weights p and q are sufficiently close to the unit for a C 1 -small f and prove our main result, Theorem 5.1, on the local uniqueness in the nonlinear problem. In the final Section 6, we investigate the question: to which extent is a symmetric tensor field f determined by data (1.10). We give a complete answer to the question in the case of M = R 3 with the Euclidean metric.

Posing the problem and introducing some notations

A smooth compact Riemannian manifold (M, g) with boundary is said to be a convex non-trapping manifold (CNTM briefly) if it satisfies two conditions: (1) the boundary ∂M is strictly convex, i.e., the second fundamental form

II(ξ, ξ) = ∇ ξ ν, ξ for ξ ∈ T x (∂M)
is positive definite for every boundary point x ∈ ∂M, where ν is the outward unit normal vector to the boundary and ∇ ξ is the covariant derivative in the direction ξ; and (2) for every x ∈ M and 0 = ξ ∈ T x M, the maximal geodesic γ x,ξ (t) determined by the initial conditions γ x,ξ (0) = x and γx,ξ (0

) = ξ is defined on a finite segment [τ -(x, ξ), τ + (x, ξ)].
In what follows, we use the notations γ x,ξ and τ ± (x, ξ) many times. They are always understood in the sense of this definition. Remark. In [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF], the term CDRM (compact dissipative Riemannian manifold) is used instead of CNTM. In the case of M ⊂ R n with the standard Euclidean metric, this definition means that M is strictly convex.

By T M = {(x, ξ) | x ∈ M, ξ ∈ T x M} we denote the tangent bundle and by

ΩM = {(x, ξ) ∈ T M | |ξ| 2 = ξ, ξ = g ij (x)ξ i ξ j = 1},
the unit sphere bundle. Its boundary can be represented as the union

∂ΩM = ∂ + ΩM ∪ ∂ -ΩM, where ∂ ± ΩM = {(x, ξ) ∈ ΩM | x ∈ ∂M, ± ξ, ν(x) ≥ 0}
is the manifold of outward (inward) unit vectors. If γ : [0, l] → M is a unit speed geodesic between boundary points, then (γ(0), γ(0)) ∈ ∂ -ΩM and (γ(l), γ(l)) ∈ ∂ + ΩM. By T C x M we denote the complexification of the tangent space T x M. The metric g determines the Hermitian scalar product on T C

x M η, ζ = g ij η i ζj .

(2.1)

For 0 = ξ ∈ T x M, by T ⊥ x,ξ M = {η ∈ T C x M | η, ξ = 0}
we denote the orthogonal complement of ξ and by π ξ : T C

x M → T C x M, the orthogonal projection onto T ⊥ x,ξ M. Let τ r s M be the bundle of complex tensors that are r times contravariant an s times covariant. Elements of the section space C ∞ (τ r s M) are smooth tensor fields of rank (r, s) on M. In the domain of a local coordinate system, such a field u ∈ C ∞ (τ r s M) can be represented by the family of smooth functions, u = (u i 1 ...ir j 1 ...js (x)), the coordinates of u, where each index takes values from 1 to n = dim M. The metric g determines canonical isomorphisms τ r s M ∼ = τ r+s 0 M ∼ = τ 0 r+s M. We will consider the isomorphisms as identifications. So, we do not distinct contra-and covariant tensors but use contra-and covariant coordinates of the same tensor. For example, for u

∈ C ∞ (τ 0 2 M) = C ∞ (τ 1 1 M) = C ∞ (τ 2 0 M), u ij = g ik u k• •j = g jk u •k i• = g ik g jl u kl .
In particular, such a tensor field determines the linear operator

u(x) : T C x M → T C x M, (uη) i = u i• •j η j
at any point x ∈ M. The product of two such operators is written in coordinates as (uv

) ij = u ik v k• •j .
The dual operator has the coordinates u * ij = ūji . The operator is Hermitian (symmetric) if and only if u ij = ūji (u ij = u ji ). The scalar product (2.1) is extended to tensors by the formula u, v = u i 1 ...im vi 1 ...ir and determines the norm

|u| 2 = u, u . For u, v ∈ C ∞ (τ 1 1 M), the norm of the product satisfies |uv| ≤ √ n|u||v|,
where n = dim M. We will also widely use semibasic tensor fields introduced in [START_REF] Pestov | Integral geometry of tensor fields on manifolds of negative curvature[END_REF], see either Section 3.4 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF] or Section 2.5 of [START_REF] Sharafutdinov | Ray transform on Riemannian manifolds (Eight lectures on integral geometry[END_REF] for a detailed presentation. Let β r s M be the bundle of complex semibasic tensor fields of rank (r, s). It is a subbundle of τ r s (T M) isomorphic to the induced bundle π * (τ r s M), where π : T M → M is the projection of the tangent bundle. A tensor u ∈ T r s,(x,ξ) (T M) at (x, ξ) ∈ T M is semibasic if it is "pure contravariant in the ξ-variable and pure covariant in x", i.e.,

u = u i 1 ...ir j 1 ...js ∂ ∂ξ i 1 ⊗ • • • ⊗ ∂ ∂ξ ir ⊗ dx j 1 ⊗ • • • ⊗ dx js .
For U ⊂ T M, by C ∞ (β r s M; U) we denote the space of smooth sections over U. The notation C ∞ (β r s M; T M) is abbreviated to C ∞ (β r s M). In the domain of a local coordinate system, such a field u ∈ C ∞ (β r s M) can be represented by the family of its coordinates, u = (u i 1 ...ir j 1 ...js (x, ξ)), which are smooth functions of 2n variables (x, ξ) = (x 1 , . . . , x n , ξ 1 , . . . , ξ n ). All the content of the previous paragraph is extended to semibasic tensor fields, where g remains the metric on M in the identification of contra-and covariant tensors. In particular, u ∈ C ∞ (β 1 1 M) determines the linear operator u(x, ξ) :

T C x M → T C x M for every (x, ξ) ∈ T M. There are two important first order differential operators v ∇ , h ∇ : C ∞ (β r s M) → C ∞ (β r s+1 M)
which are called the vertical and horizontal covariant derivatives. The operators are defined in local coordinates by the formulas v ∇k u i 1 ...ir j 1 ...js = ∂ ∂ξ k u i 1 ...ir j 1 ...js , h ∇k u i 1 ...ir j 1 ...js = ∂ ∂x k u i 1 ...ir j 1 ...js -Γ p kq ξ q ∂ ∂ξ p u i 1 .. The operator

H = ξ i h ∇i : C ∞ (β r s M) → C ∞ (β r s M
) is of the most importance in the present article. It is called the differentiation with respect to the geodesic flow.

Given a tensor field f ∈ C ∞ (τ 1 1 M) on a CNTM (M, g), let us consider the boundary value problem

HU(x, ξ) = π ξ f (x)π ξ U(x, ξ) on ΩM, U| ∂ -ΩM = E, (2.2)
where E is the identity operator. A solution U = U(x, ξ) is assumed to be a section of the bundle β 1 1 M over ΩM, i.e., U ∈ C(β 1 1 M; ΩM). In the case of M ⊂ R n with the standard Euclidean metric, f and U can be considered as n × n-matrix valued functions of x ∈ M and of (x, ξ) ∈ M × S n-1 respectively. Problem (2.2) has a unique solution. Indeed, if we restrict (2.2) to an orbit of the geodesic flow, i.e., if we set x = γ(t) and ξ = γ(t) for a unit speed geodesic γ : [0, l] → M with γ(0) ∈ ∂M, then we immediately arrive to the initial value problem (1.8). The boundary value problem (2.2) is thus equivalent to the family of initial value problems (1.8) considered for all unit speed geodesics simultaneously. The inverse problem is now formulated as follows: one has to recover the tensor field f given the trace

Φ[f ] = U| ∂ + ΩM (2.3)
of the solution to (2.2).

In order to abbreviate further formulas, let us introduce the operator P ξ on tensors which maps f (x) to π ξ f (x)π ξ for (x, ξ) ∈ ΩM, and write (2.2) in the shorter form

HU = (P ξ f )U, U| ∂ -ΩM = E. (2.4)
Because of the factor P ξ and of the boundary condition on

∂ -ΩM, the solution U to (2.4) satisfies U(x, ξ)ξ = U * (x, ξ)ξ = ξ. (2.5)
Therefore the non-trivial part of the data (2.3) consists of the restrictions Φ

[f ](x, ξ)| T ⊥ x,ξ M
for (x, ξ) ∈ ∂ + ΩM. This agrees with the above discussion of the relationship between (1.6) and (1.8). The solution U is continuous on ΩM and C ∞ -smooth on ΩM \ Ω(∂M) as one can easily prove using the strict convexity of the boundary.

Concluding the section, let us mention one more inverse problem that is not considered in the present article. Let gl(T C

x M) be the space of all linear operators on T C x M. The operator P ξ participating in (2.4) is the orthogonal projection of the space gl(T C

x M) onto the subspace

{f ∈ gl(T C x M) | f ξ = f * ξ = 0}. Let us introduce the smaller subspace {f ∈ gl(T C x M) | f ξ = f * ξ = 0, tr f = f i i = 0}
and denote by Q ξ the orthogonal projection onto the latter subspace. The corresponding inverse problem for the equation

HU = (Q ξ f )U, U| ∂ -ΩM = E (2.6)
is also of a great applied interest. To explain the physical meaning of (2.6), let us return to equation (1.5) considered in the three-dimensional case for a skew-Hermitian tensor f . The polarization vector η on (1.5) is a complex two-dimensional vector subordinate to one real condition |η| = 1. Therefore η can be described by three real parameters. Two of these parameters can be chosen to determine the shape and position of the polarization ellipse on the plane γ⊥ , while the last parameter is the phase of the electromagnetic wave. See section 6.1 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF] for a detailed discussion of the subject. Only the first two of these parameters are measured in practice. Deleting the wave phase from the data is mathematically equivalent to replacing the operator P ξ with Q ξ . The authors intend to consider the corresponding inverse problem for (2.6) in a subsequent paper.

Linear problem

Let (M, g) be a CNTM. Choose two semibasic tensor fields p, q

∈ C ∞ (β 1 1 M; ΩM) satisfy- ing p * (x, ξ)ξ = ξ, q(x, ξ)ξ = ξ. (3.1)
For a tensor field f ∈ C ∞ (τ 1 1 M), consider the boundary value problem on ΩM

Hu = p(P ξ f )q, u| ∂ -ΩM = 0. (3.2)
The problem has a unique solution u ∈ C(β 1 1 M; ΩM) and, in virtue of (3.1), the solution satisfies

u(x, ξ)ξ = u * (x, ξ)ξ = 0. (3.3)
In this section, we consider the inverse problem of recovering the tensor field f from the data

F [f ] = u| ∂ + ΩM . (3.4)
The factors p and q on (3.2) are considered as weights. We will assume the weights to be close to the unit weight E in the following sense: the inequalities

|p -E| < ε, |q -E| < ε, | v ∇ p| < ε, | v ∇ q| < ε (3.5)
hold uniformly on ΩM with the norm | • | defined in Section 2. The value of ε will be specified later. Equation (3.2) is initially considered on ΩM. To get some freedom in treating the equation, we extend it to the manifold

T 0 M = {(x, ξ) ∈ T M | ξ = 0} of nonzero vectors.
The weights are assumed to be positively homogeneous of zero degree in ξ p(x, tξ) = p(x, ξ), q(x, tξ) = q(x, ξ) for t > 0.

Then the right-hand side of (3.2) is positively homogeneous in ξ of zero degree because f is independent of ξ. The solution u must be extended to T 0 M as a homogeneous function of degree -1 u(x, tξ) = t -1 U(x, ξ) for t > 0 because the operator H increase the degree of homogeneity by 1.

Let us discuss smoothness properties of the solution u. It can be expresses by the explicit formula

u(x, ξ) = 0 τ -(x,ξ) Υ t,0 γ p(γ(t), γ(t))P γ(t) f (γ(t))q(γ(t), γ(t)) dt,
where γ = γ x,ξ and Υ t,0 γ is the parallel transport of tensors along the geodesic γ from the point γ(t) to γ(0) = x. The integrand is a smooth function. Therefore smoothness properties of u are determined by that of the integration limit τ -(x, ξ). The latter function is C ∞ -smooth on T 0 M \ T (∂M) but has singularities on T 0 (∂M). Therefore some of integrals considered below are improper and we have to verify their convergence. The verification is performed in the same way as in Section 4.6 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF]. So, in order to simplify the presentation, we will pay no attention to these singularities.

Besides (3.5), we will impose some smallness condition on the curvature of (M, g). For (x, ξ) ∈ ΩM, let K(x, ξ) be the supremum of the absolute values of sectional curvatures at the point x over all two-dimensional subspaces of and for any weights p, q ∈ C ∞ (β 1 1 M; ΩM) satisfying (3.1) and (3.5), every tensor field f ∈ C ∞ (τ 1 1 M) can be uniquely recovered from the trace (3.4) of the solution to the boundary value problem (3.2) and the stability estimate

T x M containing ξ. Define k(M, g) = sup (x,ξ)∈∂ -ΩM τ + (x,ξ) 0 tK(γ x,ξ (t), γx,ξ (t))dt. ( 3
f L 2 ≤ C F [f ] H 1 (3.8)
holds with a constant C independent of f . In the case of n = 3, the same statement is true for a symmetric tensor field f .

In the case of a real symmetric f and unit weights, this theorem is a partial case of Theorem 5.2.2 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF]. We will show that the same proof works with some modifications for Theorem 3.1.

Proof of Theorem 3.1. We rewrite equation (3.2) in the form

Hu = P ξ f + r, (3.9)
where r = (p -E)P ξ f + pP ξ f (q -E).

(3.10)

The remainder r is small by (3.5). Because of (3.3), the function

u = (u ij (x, ξ)) is orthogonal to ξ in both indices ξ i u ij = ξ j u ij = 0. (3.11)
We write down the Pestov identity for the semibasic tensor field u (see Lemma 4.4.1 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF] for the case of a real u and Lemma 5.1 of [START_REF] Sharafutdinov | Inverse problem of determining a connection on a vector bundle[END_REF] for the general case)

2 Re h ∇ u, v ∇ Hu = | h ∇ u| 2 + h ∇i v i + v ∇i w i -R 1 [u],
(3.12) where •, • and | • | are the scalar product and norm on semibasic tensors defined in Section 2,

v i = Re ξ i h ∇ j u i 1 i 2 • v ∇j ūi 1 i 2 -ξ j v ∇ i u i 1 i 2 • h ∇j ūi 1 i 2 , (3.13)
w i = Re ξ j h ∇ i u i 1 i 2 • h ∇j ūi 1 i 2 , (3.14) R 1 [u] = R kplq ξ p ξ q v ∇ k u i 1 i 2 • v ∇ l ūi 1 i 2 + Re (R i 1 pqj u pi 2 + R i 2 pqj u i 1 p )ξ q v ∇ j ūi 1 i 2 , (3.15)
and (R ijkl ) is the curvature tensor. Identity (3.12) holds on any open D ⊂ T M for every u ∈ C 2 (β 1 1 M; D). In our case, D = T 0 M \ T (∂M) and u ∈ C ∞ (β 1 1 M; D) as is shown above.

The most part of the proof deals with the left-hand side of (3.12). We will first transform it by distinguishing some divergent terms and then will estimate it.

From (3.9)

h ∇ u, v ∇ Hu = h ∇ u, v ∇ (P ξ f ) + h ∇ u, v ∇ r . (3.16)
We will first investigate the first term on the right-hand side of (3.16). To this end we represent f as

f ij (x) = fij (x, ξ) + ξ j a i (x, ξ) + ξ i bj (x, ξ) + ξ i ξ j c(x, ξ), (3.17)
where ( fij ) is a semibasic tensor field orthogonal to ξ in both indices

fij ξ i = fij ξ j = 0, (3.18)
semibasic covector fields a and b are orthogonal to ξ

a i ξ i = b i ξ i = 0, (3.19)
and c(x, ξ) is a scalar function. One can easily check the existence and uniqueness of the representation. The (vector versions of the) fields a and b are expressed through f by the formulas

a = 1 |ξ| 2 π ξ f ξ, b = 1 |ξ| 2 π ξ f * ξ.
(3.20)

As follows from (3.17)-(3.19),

P ξ f = f or in coordinates (P ξ f ) ij = fij = f ij -a i ξ j -bj ξ i -cξ i ξ j . (3.21)
Differentiating the last equality with respect to ξ and using the fact that f is independent of ξ, we obtain

v ∇k (P ξ f ) ij = -ξ j v ∇k a i -ξ i v ∇k bj -ξ i ξ j v ∇k c -g jk a i -g ik bj -(g ik ξ j + g jk ξ i )c. Therefore h ∇ u, v ∇ (P ξ f ) = h ∇ k u ij • v ∇k (P ξ f) ij = = h ∇ k u ij -ξ j v ∇k āi -ξ i v ∇k b j -ξ i ξ j v ∇k c -g jk āi -g ik b j -(g ik ξ j + g jk ξ i )c .
The tensor h ∇ k u ij is orthogonal to ξ in the indices i and j as follows from (3.11). Therefore the last formula is simplified to the following one:

h ∇ u, v ∇ (P ξ f ) = - h ∇ p u ip • āi - h ∇ p u pi • b i .
Introducing the semibasic covector fields h δ1 u and h δ2 u by the equalities

( h δ1 u) i = h ∇ p u ip , ( h δ2 u) i = h ∇ p u pi , (3.22)
we write the result in the form

h ∇ u, v ∇ (P ξ f ) = - h δ1 u, a - h δ2 u, b . (3.23)
This implies the estimate

2 Re h ∇ u, v ∇ (P ξ f ) ≤ β 2 (| h δ1 u| 2 + | h δ2 u| 2 ) + 2 β (|a| 2 + |b| 2 ), (3.24)
where β is an arbitrary positive number.

Next, we transform the expression | h δ1 u| 2 by distinguishing a divergent term

| h δ1 u| 2 = ( h δ1 u) i ( h δ1 ū) i = h ∇p u ip • h ∇ q ūiq = = h ∇p (u ip h ∇ q ūiq ) -u ip h ∇p h ∇ q ūiq = h ∇p (u ip h ∇ q ūiq ) -u •p i• h ∇p h ∇q ūiq . (3.25)
By the commutator formula for horizontal derivatives (see Theorem 3.5.2 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF]),

h ∇p h ∇q ūiq = h ∇q h ∇p ūiq -R k jpq ξ j v ∇k ūiq + R i jpq ūjq + R q jpq ūij .
Substituting this value into the previous formula, we obtain

| h δ1 u| 2 = Re -u •p i• h ∇q h ∇p ūiq + h ∇p (u ip h ∇ q ūiq ) + R 2 [u], (3.26)
where

R 2 [u] = Re u •p i• (R k jpq ξ j v ∇k ūiq -R i jpq ūjq -R q jpq ūij ) .
We now transform the first summand on the right-hand side of (3.26) in the order reverse to that used in (3.25)

| h δ1 u| 2 = Re - h ∇q (u •p i• h ∇p ūiq ) + h ∇q u •p i• • h ∇p ūiq + h ∇p (u ip h ∇ q ūiq ) + R 2 [u].
Introducing the semibasic vector field ṽ1 by the formula

(ṽ 1 ) i = Re u ji h ∇ k ūjk -u jk h ∇ k ūji , (3.27)
we write the result in the form

| h δ1 u| 2 = Re h ∇ i u jk • h ∇k ūji + h ∇i (ṽ 1 ) i + R 2 [u]. (3.28)
In the same way, we obtain

| h δ2 u| 2 = Re h ∇ i u jk • h ∇j ūik + h ∇i (ṽ 2 ) i + R 3 [u] (3.29) with (ṽ 2 ) i = Re u ij h ∇ k ūkj -u kj h ∇ k ūij (3.30) and R 3 [u] = Re u p• •i (R k jpq ξ j v ∇k ūqi -R q jpq ūji -R i jpq ūqi ) .
Taking the sum of (3.28) and (3.29), we have

| h δ1 u| 2 + | h δ2 u| 2 = Re h ∇ i u jk • h ∇j ūik + h ∇ i u jk • h ∇k ūji + h ∇i ṽi + R 4 [u], (3.31)
where ṽ

= ṽ1 + ṽ2 , R 4 [u] = R 2 [u] + R 3 [u].
Introduce the semibasic tensor field z = (z ijk ) by the formula

h ∇i u jk = ξ i |ξ| 2 (Hu) jk + z ijk .
(3.32)

The idea of this new notation is that the tensor z is orthogonal to ξ in all its indices

ξ i z ijk = ξ j z ijk = ξ k z ijk = 0, ( 3
| h ∇ u| 2 = 1 |ξ| 2 |Hu| 2 + |z| 2 .
(3.34)

The first two terms on the right-hand side of (3.31) can be expressed through z. Indeed, one easily see with the help of (3.32) and (3.33) that Re

h ∇ i u jk • h ∇j ūik + h ∇ i u jk • h ∇k ūji = Re (z ijk zjik + z ijk zkji ) ≤ 2|z| 2 .
With the help of the last inequality, (3.31) implies the estimate

| h δ1 u| 2 + | h δ2 u| 2 ≤ 2|z| 2 + h ∇i ṽi + R 4 [u] (3.35)
which, together with (3.24), gives

2 Re h ∇ u, v ∇ (P ξ f ) ≤ β|z| 2 + 2 β (|a| 2 + |b| 2 ) + β 2 h ∇i ṽi + β 2 R 4 [u].
Substitute values (3.20) for a and b into the last formula

2 Re h ∇ u, v ∇ (P ξ f ) ≤ β|z| 2 + 2 β|ξ| 4 (|π ξ f ξ| 2 + |π ξ f * ξ| 2 ) + β 2 h ∇i ṽi + β 2 R 4 [u]. (3.36)
Next, we estimate the second term on the right-hand side of (3.16). We differentiate equality (3.10) with respect to ξ taking the independence f of ξ into account

v ∇ r = v ∇ ((p -E)P ξ ) f + v ∇ (pP ξ )f (q -E) + pP ξ f v ∇ q.
In what follows in the proof, we denote different constants depending only on n = dim M by the same letter C. On using (3.5), we obtain from the last formula

| v ∇ r| ≤ Cε |ξ| |f |. From this 2 Re h ∇ u, v ∇ r ≤ Cε(| h ∇ u| 2 + 1 |ξ| 2 |f | 2 ).
(3.37)

Combining (3.36) and (3.37), we obtain from (3.16)

2 Re h ∇ u, v ∇ Hu ≤ β|z| 2 + 2 β|ξ| 4 (|π ξ f ξ| 2 +|π ξ f * ξ| 2 )+ β 2 h ∇i ṽi +Cε(| h ∇ u| 2 + 1 |ξ| 2 |f | 2 )+ β 2 R 4 [u].
(3.38) Estimating the left-hand side of the Pestov identity (3.12) by (3.38), we obtain for |ξ| = 1

| h ∇ u| 2 + v ∇i w i -β|z| 2 - 2 β (|π ξ f ξ| 2 +|π ξ f * ξ| 2 )-Cε(| h ∇ u| 2 +|f | 2 ) ≤ h ∇i ( β 2 ṽi -v i )+R[u], (3.39)
where

R[u] = R 1 [u] + β 2 R 4 [u].
We multiply inequality (3.39) by the volume form dΣ = |d x ω(ξ) ∧ dV n (x)|, integrate over ΩM, and transform the integrals of divergent terms by Gauss-Ostrogradskii formulas (see Theorem 3.6.3 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF])

ΩM | h ∇ u| 2 + (n -2)|Hu| 2 -β|z| 2 - 2 β (|π ξ f ξ| 2 + |π ξ f * ξ| 2 ) -Cε(| h ∇ u| 2 + |f | 2 ) dΣ ≤ ≤ ∂ΩM β 2 ṽ -v, ν dΣ 2n-2 + ΩM R[u]dΣ. (3.40)
The second term on the left-hand side has appeared because w is positively homogeneous of degree -1 in ξ and satisfies ξ, w = |Hu| 2 as is seen from (3.14). Substituting the value | h ∇ u| 2 = |z| 2 + |Hu| 2 from (3.34), we write the result in the form

ΩM (1 -β -Cε)|z| 2 + (n -1 -Cε)|Hu| 2 - 2 β (|π ξ f ξ| 2 + |π ξ f * ξ| 2 ) -Cε|f | 2 dΣ ≤ ≤ ∂ΩM β 2 ṽ -v, ν dΣ 2n-2 + ΩM R[u]dΣ. (3.41)
Assuming β ≤ 1, integrals on the right-hand side of (3.41) can be estimated exactly as in Section 5.5 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF] ΩM

R[u]dΣ ≤ Ck(M, g) ΩM | h ∇ u| 2 dΣ ≤ Cδ ΩM | h ∇ u| 2 dΣ, (3.42) ∂ΩM β 2 ṽ -v, ν dΣ 2n-2 ≤ D u| ∂ + ΩM 2 H 1 , (3.43)
where k(M, g) is defined by (3.6). The second inequality on (3.42) is valid because of (3.7). The constant D on (3.43) depends on (M, g), unlike the constant C on (3.42) which depends only on n.

Combining (3.41)-(3.43) and using again the equality

| h ∇ u| 2 = |z| 2 + |Hu| 2 , we obtain ΩM (1 -β -Cε -Cδ)|z| 2 + (n -1 -Cε -Cδ)|Hu| 2 - 2 β (|π ξ f ξ| 2 + |π ξ f * ξ| 2 ) -Cε|f | 2 dΣ ≤ ≤ D u| ∂ + ΩM 2 H 1 (3.44)
with some new constant C depending only on n.

Let us compare |Hu| and |P ξ f |. The estimate |r| < Cε|f | follows from (3.5) and (3.10). The latter, together with (3.9), implies

|Hu| 2 ≥ |P ξ f | 2 -Cε|f | 2 .
(3.45)

Using the last inequality, we transform (3.44) to the final form

M ΩxM (1 -β -Cε -Cδ)|z| 2 + (n -1 -Cε -Cδ)|P ξ f | 2 - - 2 β (|π ξ f ξ| 2 + |π ξ f * ξ| 2 ) -Cε|f | 2 dω x (ξ)dV n (x) ≤ D u| ∂ + ΩM 2 H 1 . (3.46)
Let us remind that β is an arbitrary number satisfying 0 < β ≤ 1.

Lemma 3.2 For every Riemannian manifold (M, g) of dimension n ≥ 4 and every point x ∈ M, the Hermitian form

B(f, f ) = ΩxM (n -1)|P ξ f | 2 -2(|π ξ f ξ| 2 + |π ξ f * ξ| 2 ) dω x (ξ)
is positive definite on the space of second rank tensors at x. Moreover, the estimate

B(f, f ) ≥ c|f | 2
holds with a positive constant c depending only on n. In the case of n = 3, the same statement is valid for symmetric tensors.

The proof of the lemma will be given later, and now we finish the proof of Theorem 3.1 by making use of the lemma.

By the lemma, the inequality

c f 2 L 2 = c M |f | 2 dV n (x) ≤ M ΩxM (1 -β -Cε -Cδ)|z| 2 + +(n -1 -Cε -Cδ)|P ξ f | 2 - 2 β (|π ξ f ξ| 2 + |π ξ f * ξ| 2 ) -Cε|f | 2 dω x (ξ)dV n (x) (3
f 2 L 2 ≤ C u| ∂ + ΩM 2 H 1 = C F [f ] 2 H 1
with C = D/c. This finishes the proof of Theorem 3.1.

Proof of Lemma 3.2. One easily checks the equality B(f, f ) = B(u, u) + B(v, v) for a complex tensor f represented as f = u + iv with real u and v. Therefore it suffices to prove the statement for a real tensor f .

The corresponding symmetric bilinear form is

B(f, h) = ΩxM (n -1) P ξ f, h -2( π ξ f ξ, hξ + π ξ f * ξ, h * ξ ) dω x (ξ).
Obviously, B(f, h) = 0 for a symmetric f and skew-symmetric h. Therefore, it suffices to prove the positiveness of B on the spaces of real symmetric and skew-symmetric tensors separately.

The positiveness of B on the space of real symmetric tensors is proved in Lemma 5.6.1 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF], where π ξ f ξ is denoted by P ξ j ξ f . Thus, we have to consider the quadratic form B on the space of real skew-symmetric tensors. On making use of an orthonormal basis, we identify T x M with R n endowed with the standard scalar product and identify Ω x M with the unit sphere Ω ⊂ R n . So

1 ω B(f, f ) = 1 ω Ω (n -1)|P ξ f | 2 -4|π ξ f ξ| 2 dω(ξ), (3.48)
where ω is the volume of Ω and dω is the standard volume form on Ω. For a skew-symmetric f π ξ f ξ = f ξ (3.49) since f ξ is orthogonal to ξ. Therefore formulas (3.20) and (3.21) are simplified to the following ones:

a = -b = f ξ, (P ξ f ) ij = f ij -(f ξ) i ξ j + ξ i (f ξ) j for |ξ| = 1.
Using the last equality and f ξ, ξ = 0, we easily calculate

|P ξ f | 2 = |f | 2 -2|f ξ| 2 for |ξ| = 1. ( 3
1 ω B(f, f ) = (n -1)|f | 2 - 2(n + 1) ω Ω |f ξ| 2 dω. (3.51)
On using the obvious relation

1 ω Ω ξ i ξ j dω = 0 for i = j, 1/n for i = j, we find 1 ω Ω |f ξ| 2 dω = 1 n |f | 2 .
Inserting this value into (3.51), we see that

1 ω B(f, f ) = n 2 -3n -2 n |f | 2 .
This implies the positiveness of B on skew-symmetric tensors for n ≥ 4.

Three-dimensional case

We will first show that both our problems, linear and nonlinear, possess some nonuniqueness in the three-dimensional case. Let (M, g) be a three-dimensional CNTM which is assumed to be oriented. Every tangent space T x M is a three-dimensional oriented Euclidean space. So, the vector product

T x M × T x M → T x M, (v, w) → v × w is well defined. It is extended to the C-bilinear operation T C x M × T C x M → T C x M, (v, w) → v × w on complex vectors. For a complex vector field v ∈ C ∞ (τ 1 0 M), we denote by L v ∈ C ∞ (τ 1 1 M) the operator of vector multiplication by v, L v (x)η = v(x) × η for η ∈ T C x M.
Note that L v is a skew-symmetric tensor field. Quite similarly, for a semibasic vector field v ∈ C ∞ (β 1 0 M), the operator L v ∈ C ∞ (β 1 1 M) is defined. Let us prove the formula

P ξ L v = π ξ L v π ξ = v, ξ |ξ| 2 L ξ . (4.1)
We remind that •, • and | • | are defined in Section 2. This is a pure algebraic local formula. So, we can use a positive orthonormal basis (e 1 , e 2 , e 3 = ξ/|ξ|) in T x M. In such a basis

π ξ =   1 0 0 0 1 0 0 0 0   , L v =   0 -v 3 v 2 v 3 0 -v 1 -v 2 v 1 0   , L ξ =   0 -|ξ| 0 |ξ| 0 0 0 0 0   , v, ξ = v 3 |ξ|,
and the formula follows immediately. Next, we prove the formula

h ∇ L ξ = 0. (4.2)
The formula is quite expectable since h ∇ ξ = 0. Nevertheless, it needs a proof. Here, we have to use a general coordinate system since h ∇ is a differential operator. The vector product is expressed by the formula

(v × w) i = 1 √ g (v i+1 w i+2 -v i+2 w i+1 )
in general coordinates, where g = det(g ij ) and indices are reduced modulo 3. Therefore

(L v ) i,i+1 = - √ gv i+2 (4.3)
and (4.2) is proved by straightforward calculations with making use of ∇ i g jk = 0 and h ∇i ξ j = 0. We derive from (4.1)-(4.2) and the formula H(1/|ξ| 2 ) = 0 that For a real function λ, U(x, ξ) is the rotation of T x M around the axis ξ by the angle λ(x).

H( λ |ξ| 2 L ξ ) = Hλ |ξ| 2 L ξ = P ξ L ∇λ for a complex function λ ∈ C ∞ (M). Thus H( λ |ξ| 2 L ξ ) = P ξ L ∇λ . ( 4
In the case of a complex λ, the operator is also well defined although its geometric sense is more complicated. The semibasic tensor field U ∈ C ∞ (β where λ(t) = λ(γ(t)) and U(t) = U(γ(t), γ(t)). Conversely, if (4.10) holds for any unit speed geodesic γ, then (4.9) is true. To prove (4.10), we choose an orthonormal basis (e 1 (t), e 2 (t), e 3 (t) = γ(t)) of T γ(t) M which is parallel along γ. In such a basis, (4.10) is equivalent to the matrix equation

d dt   cos λ(t) -sin λ(t) 0 sin λ(t) cos λ(t) 0 0 0 1   = dλ dt   0 -1 0 1 0 0 0 0 0     cos λ(t) -sin λ(t) 0 sin λ(t) cos λ(t) 0 0 0 1  
which is obviously true. Comparing (4.7)-(4.8) with (2.3)-(2.4), we see that a field f = L ∇λ with λ| ∂M = 0 cannot be recovered from data (2.3).

Let us introduce the definition:

f ∈ C ∞ (τ 1 1 M
) is said to be a potential field if it can be represented as f = L ∇λ for some function λ ∈ C ∞ (M) vanishing on the boundary, λ| ∂M = 0. Potential fields constitute the natural obstruction for the uniqueness in the both, linear and nonlinear, problems. We are going to prove that a solution to the linear problem is unique up to the obstruction. The corresponding local result for the nonlinear problem will be obtained in the next section. First of all we prove Lemma 4.1 (on decomposition). Let (M, g) be a compact oriented three-dimensional Riemannian manifold with boundary. Every tensor field f ∈ C ∞ (τ 1 1 M) can be uniquely represented as

f = L ∇λ + f (4.11)
with some λ ∈ C ∞ (M) satisfying λ| ∂M = 0 (4.12)

and some tensor field

f ∈ C ∞ (τ 1 1 M) satisfying the condition: the 2-form fij dx i ∧ dx j is closed, d( fij dx i ∧ dx j ) = 0. (4.13)
The summands of decomposition (4.11) are called the potential and closed parts of f respectively. Note that (4.13) involves only the skew-symmetric part of f , i.e., a symmetric tensor field is closed. Lemma 4.1 can be derived from the Hodge-Morrey decomposition [START_REF] Schwarz | Hodge Decomposition -a Method for Solving Boundary Value Problems[END_REF] but we give a shorter independent proof.

Proof of Lemma 4.1. We first prove the uniqueness statement. Assume (4.11)-(4.13) to be valid. Applying the exterior derivative d to the form f ij dx i ∧ dx j and using (4.11) and (4.13), we obtain

d (L ∇λ ) ij dx i ∧ dx j = d(f ij dx i ∧ dx j ).
On using (4.3), one can check by a straightforward calculation in coordinates that

d (L ∇λ ) ij dx i ∧ dx j = -2∆λ √ gdx 1 ∧ dx 2 ∧ dx 3 ,
where ∆ is the Laplace-Beltrami operator. Thus, the function λ solves the Dirichlet problem

∆λ = - 1 √ g ∂f - 23 ∂x 1 + ∂f - 31 ∂x 2 + ∂f - 12 ∂x 3 , λ| ∂M = 0, (4.14)
where f -is the skew-symmetric part of f , i.e., f - ij = 1 2 (f ij -f ji ). The solution to the Dirichlet problem is unique.

The existence statement is proved by reverse arguments. Given f , let λ be the solution to the Dirichlet problem (4.14) and f = f -L ∇λ . Then (4.11)-(4.13) holds. The lemma is proved.

We restrict ourselves to considering the inverse problem for closed fields only.

Theorem 4.2 There exist such positive numbers δ and ε that, for any oriented 3-dimensional CNTM (M, g) satisfying (3.7) and for any weights p, q ∈ C ∞ (β 1 1 M; ΩM) satisfying (3.1) and (3.5), every closed tensor field f ∈ C ∞ (τ 1 1 M) can be uniquely recovered from the trace (3.4) of the solution to the boundary value problem (3.2) and the stability estimate

f 2 L 2 ≤ C F [f ] 2 H 1 + f | ∂M L 2 • F [f ] L 2 (4.15)
holds with a constant C independent of f .

Proof follows the same line as the proof of Theorem 3.1 with the following difference: the left-hand side of the Pestov identity (3.12) will be estimated in a different way by making use of the closeness of f . We represent f as the sum of symmetric and skew-symmetric fields

f = f + + f -, f + ij = f + ji , f - ij = -f - ji .
Taking the symmetry of Christoffel symbols into account, the closeness condition for f can be written as

∇ i f - jk + ∇ j f - ki + ∇ k f - ij = 0. (4.16)
The vector f -ξ is orthogonal to ξ and therefore π ξ f -ξ = f -ξ. Formulas (3.20) take now the form

a = 1 |ξ| 2 π ξ f + ξ + 1 |ξ| 2 f -ξ, b = 1 |ξ| 2 π ξ f + ξ - 1 |ξ| 2 f -ξ. (4.17)
We write equation (3.2) in the form (3.9) with the remainder r defined by (3.10). Then we write the Pestov identity (3.12) for u with terms defined by (3.13)-(3.15). The left hand-side of the identity can be written as in (3.16).

The main problem is estimating the first term h ∇ u, v ∇ (P ξ f ) on the right-hand side of (3.16). To this end we represent f in form (3.17) with a and b defined by (4.17). Then (3.23) holds. In view of (4.17), equation (3.23) can be written as

h ∇ u, v ∇ (P ξ f ) = - 1 |ξ| 2 h δ1 u + h δ2 u, π ξ f + ξ - 1 |ξ| 2 h δ1 u - h δ2 u, f -ξ . (4.18)
We transform the second term on the right-hand side of (4.18) by distinguishing divergent terms. Using definition (3.22) of h δ1 u and h δ2 u, we write

h δ1 u - h δ2 u, f -ξ = ( h ∇p u ip - h ∇p u pi ) f - ik ξ k = = h ∇p (ξ k u ip f- ik + ξ k u pi f- ki ) -ξ k u ip ( h ∇p f - ik + h ∇i f - kp ).
Now, we use the closeness condition (4.16) to transform this formula to the following one:

h δ1 u - h δ2 u, f -ξ = ξ k u ip h ∇k f - pi + h ∇p (ξ k u ip f - ik + ξ k u pi f - ki )
. Finally, we distinguish a divergent term from the first summand on the right-hand side

h δ1 u - h δ2 u, f -ξ = h ∇k (ξ k u ip f - pi ) -ξ k h ∇k u ip • f - pi + h ∇p (ξ k u ip f - ik + ξ k u pi f- ki )
. This can be written as

h δ1 u - h δ2 u, f -ξ = Hu, f -+ h ∇i (ξ i u jk f - kj + ξ k u ji f - jk + ξ k u ij f - kj ). (4.19)
Next, we calculate the first term on the right-hand side of (4.19) by making use of (3.9) and (3.21)

Hu, f -= (Hu) ij f- ij = ((P ξ f ) ij + r ij ) f - ij = (f ij -a i ξ j -bj ξ i -cξ i ξ j + r ij ) f - ij . (4.20)
Since f -is a skew-symmetric tensor, cξ i ξ j fij = 0. This means that the term cξ i ξ j on the right-hand side of (4.20) can be omitted. By the same reason, the term f ij can be replaced with (f -) ij , the vector a = 1 |ξ| 2 π ξ f ξ can be replaced with 1 |ξ| 2 f ξ, and the vector b can be replaced with 1 |ξ| 2 f * ξ. In such the way, (4.20) takes the form

Hu, f -= |f -| 2 - 1 |ξ| 2 (f ik ξ k ξ j + f kj ξ k ξ i ) f - ij + r, f -.
The second term on the right-hand side is independent of f + , and the formula can be written as

Hu, f -= |f -| 2 - 2 |ξ| 2 |f -ξ| 2 + r, f -. (4.21)
Substitute (4.21) into (4.19)

h δ1 u - h δ2 u, f -ξ = |f -| 2 - 2 |ξ| 2 |f -ξ| 2 + r, f -+ h ∇i (ξ i u jk f- kj + ξ k u ji f- jk + ξ k u ij f- kj ).
Inserting this expression into (4.18) and estimating the first term on the right-hand side of (4.18) in a similar way as in deriving (3.24), we arrive to the inequality

2 Re h ∇ u, v ∇ (P ξ f ) ≤ β 2 (| h δ1 u| 2 + | h δ2 u| 2 ) + 4 β|ξ| 4 |π ξ f + ξ| 2 - - 2 |ξ| 2 | f-| 2 + 4 |ξ| 4 |f -ξ| 2 + h ∇i (ṽ 3 ) i + 2 |ξ| 2 |r| • |f |, (4.22)
where (ṽ 3

) i = - 2 |ξ| 2 Re(ξ i u jk f - kj + ξ k u ji f - jk + ξ k u ij f - kj ) (4.
∇ u, v ∇ (P ξ f ) ≤ β|z| 2 + 4 β|ξ| 4 |π ξ f + ξ| 2 - 2 |ξ| 2 | f-| 2 + 4 |ξ| 4 |f -ξ| 2 + Cε |ξ| 2 |f | 2 + h ∇i ṽi + β 2 R 4 [u] ( 4
∇ u, v ∇ Hu ≤ β|z| 2 + 4 β |π ξ f + ξ| 2 -2| f -| 2 + 4|f -ξ| 2 + Cε(|f | 2 + | h ∇ u| 2 ) + h ∇i ṽi + β 2 R 4 [u] (4.25) for |ξ| = 1.
We estimate the left-hand side of the Pestov identity (3.12) by (4.25) and write the result in the form

(1-Cε)| h ∇ u| 2 + v ∇i w i -β|z| 2 - 4 β |π ξ f + ξ| 2 +2| f -| 2 -4|f -ξ| 2 -Cε|f | 2 ≤ h ∇i (ṽ i -v i )+ β 2 R 4 [u].
Integrating this inequality and transforming the integrals of divergent terms in the same way as in (3.40), we obtain

ΩM (1 -Cε)| h ∇ u| 2 + |Hu| 2 -β|z| 2 - 4 β |π ξ f + ξ| 2 + 2| f -| 2 -4|f -ξ| 2 -Cε|f | 2 dΣ ≤ ≤ ∂ΩM ṽ -v, ν dΣ 2n-2 + ΩM R 4 [u]dΣ.
Substituting the value | h ∇ u| 2 = |z| 2 + |Hu| 2 from (3.34), we write the result in the form

ΩM (1 -β -Cε)|z| 2 + (2 -Cε)|Hu| 2 - 4 β |π ξ f + ξ| 2 + 2| f-| 2 -4|f -ξ| 2 -Cε|f | 2 dΣ ≤ ≤ ∂ΩM ṽ -v, ν dΣ 2n-2 + ΩM R 4 [u]dΣ. (4.26)
The curvature dependent integral on (4.26) is estimated as before in (3.42

) ΩM R 1 [u]dΣ ≤ Cδ ΩM | h ∇ u| 2 dΣ, (4.27)
while the boundary integral is estimated in a little bit different way. Namely, instead of (3.43), we have the estimate

∂ΩM ṽ -v, ν dΣ 2n-2 ≤ D u| ∂ + ΩM 2 H 1 + f | ∂M L 2 • u| ∂ + ΩM L 2 . (4.28)
The second term on the right-hand side of (4.28) appears because of the dependence of ṽ on f as is seen from (4.23). Inequality (4.28) is proved in the same way as estimate (4.7.2) of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF].

Combining (4.26)-(4.28) and using again the equality

| h ∇ u| 2 = |z| 2 + |Hu| 2 , we obtain ΩM (1-β-Cδ-Cε)|z| 2 +(2-Cδ-Cε)|Hu| 2 - 4 β |π ξ f + ξ| 2 +2| f-| 2 -4|f -ξ| 2 -Cε|f | 2 dΣ ≤ ≤ D u| ∂ + ΩM 2 H 1 + f | ∂M L 2 • u| ∂ + ΩM L 2 . (4.29)
The tensors f + and f -are orthogonal to each other as well as P ξ f + and P ξ f -are orthogonal to each other. Therefore

|f | 2 = |f + | 2 + |f -| 2 , |P ξ f | 2 = |P ξ f + | 2 + |P ξ f -| 2 .
With the help of these equalities, (3.45) gives

|Hu| 2 ≥ |P ξ f + | 2 -Cε|f + | 2 -Cε|f -| 2 .
(4.30)

We use (4.30) to transform the estimate (4.29) to the final form

(1 -β -Cδ -Cε) ΩM |z| 2 dΣ + ΩM (2 -Cδ -Cε)|P ξ f + | 2 - 4 β |π ξ f + ξ| 2 -Cε|f + | 2 dΣ+ +2 ΩM (| f -| 2 -2|f -ξ| 2 -Cε|f -| 2 )dΣ ≤ D u| ∂ + ΩM 2 H 1 + f | ∂M L 2 • u| ∂ + ΩM L 2 (4.31)
with some new constant C. For β = 1 and ε = δ = 0, the left-hand side of (4.31) is the integral over M of the Hermitian form

A(f, f ) = B(f + , f + ) + 2Q(f -, f -),
where B is the same as in Lemma 3.2 and

Q(f, f ) = ΩxM (|f | 2 -2|f ξ| 2 )dω x (ξ). (4.32)
By Lemma 3.2, the form B is positively definite on symmetric tensors in the 3-dimensional case. The form Q is positively definite on skew-symmetric tensors. Indeed, repeating the arguments of the proof of Lemma 3.2, we derive in the 3-dimensional case

1 ω ΩxM (|f | 2 -2|f ξ| 2 )dω x (ξ) = 1 3 |f | 2
for a skew-symmetric tensor f . Now, the proof is finished in the same way as in Theorem 3.1.

Nonlinear problem

We return to considering the inverse problem of recovering a tensor field f ∈ C ∞ (τ 1 1 M) on a CNTM (M, g) from the data

Φ[f ] = U| ∂ + ΩM , (5.1)
where U ∈ C(β 1 1 M; ΩM) is the solution to the boundary value problem on ΩM

HU = (P ξ f )U, U| ∂ -ΩM = E.
(5.2)

We will prove the uniqueness under the following smallness assumptions on f :

|f (x)| < ε for x ∈ ∂M, (5.3
)

0 τ -(x,ξ) |f (γ x,ξ (t))|dt < ε, 0 τ -(x,ξ) |∇f (γ x,ξ (t))|dt < ε for (x, ξ) ∈ ∂ + ΩM. ( 5.4)
We remind that we use notations introduced in Section 2. In particular, ∇ is the covariant derivative. Note that these smallness conditions are quite similar to that of Theorem 2 of [START_REF] Vertgeim | Integral geometry with matrix weight, and a nonlinear problem of recovering matrices[END_REF].

Theorem 5.1 It is possible to choose a positive number δ = δ(n) for n ≥ 4 such that, for an n-dimensional CNTM (M, g) satisfying the curvature condition (3.7), there exists a positive number ε = ε(M, g) such that the following statement is true. Let two tensor fields

f i ∈ C ∞ (τ 1 1 M) (i = 1, 2) satisfy (5.
3)-(5.4) and Φ i = Φ[f i ] be the corresponding data. Then the estimate

f 2 -f 1 L 2 ≤ C Φ -1 1 Φ 2 -E H 1 (5.5)
holds with a constant C independent of f i . In particular,

f 1 = f 2 if Φ 1 = Φ 2 .
In the case of n = 3, the same statement is true under the additional assumption that f 2 -f 1 is a closed tensor field.

Proof. Let U i ∈ C(β 1 1 M; ΩM) (i = 1, 2) be the solution to the boundary value problem

HU i = (P ξ f i )U i , U i | ∂ -ΩM = E. (5.6)
According to (5.2), the solution satisfies

U i ξ = U * i ξ = ξ. (5.7) Set u = U -1 1 U 2 -E. Using the equalities HE = 0, H(U -1 1 U 2 ) = HU -1 1 • U 2 + U -1 1 HU 2 , and HU -1 1 = -U -1 1 HU 1 • U -1
1 , one easily derive from (5.6) that u solves the boundary value problem

Hu = U -1 1 P ξ (f 2 -f 1 ) U 2 , u| ∂ -ΩM = 0 (5.8) and u| ∂ + ΩM = Φ -1 1 Φ 2 -E.
(5.9)

Setting f = f 2 -f 1 , p = U - 1
1 , q = U 2 , we write (5.8)-(5.9) in the form Hu = p(P ξ f )q, u| ∂ -ΩM = 0, (5.10)

F [f ] = u| ∂ + ΩM = Φ -1 1 Φ 2 -E.
(5.11)

We have arrived to the linear problem considered in Sections 3 and 4. If we will prove that the weights p = U -1 1 and q = U 2 satisfy conditions (3.1) and (3.5), we would be able to apply Theorems 3. with some constant C depending on (M, g) but not on f .

To prove Lemma 5.2, we need the following estimate for solutions to linear ordinary differential equations (see Lemma 4.1 of Chapter IV of [START_REF] Ph | Ordinary Differential Equations[END_REF]): The constant C appears in Lemma 5.4 since different norms are used in this Lemma and Lemma 5. [START_REF] Yu | Geometric Optics of Inhomogeneous Media[END_REF]. In what follows in this section, we denote different constants depending on (M, g) by the same letter C.

M; ∂ -ΩM ). Then |y(x, ξ)| ≤ C |y 0 (γ(τ -(x, ξ)), γ(τ -(x, ξ)))| + 0 τ -(x,ξ) |f (γ(t), γ(t))|dt × × exp C 0 τ -(x,ξ) |A(γ(t), γ (
Proof of Lemma 5.2. Let U be the solution to (5.2). Then U -E solves the boundary value problem

H(U -E) = (P ξ f )(U -E) + P ξ f, (U -E)| ∂ -ΩM = 0.
Applying Lemma 5.4, we obtain the estimate

|(U -E)(x, ξ)| ≤ C 0 τ -(x,ξ) |P γ(t) f (γ(t))|dt exp C 0 τ -(x,ξ) |P γ(t) f (γ(t))|dt ≤ ≤ C 0 τ -(x,ξ) |f (γ(t))|dt exp C 0 τ -(x,ξ) |f (γ(t))|dt ,
Together with (5.4), this implies the first of inequalities (5.12).

The proof of the second of estimates (5.12) is more troublesome because

U i j )| x n =0 = (P ξ f ) i j . Since ξ n = -ξ, ν
, where ν = ν(x) is the unit outward normal to the boundary, this gives with the help of (5.3) ∇ and H is a little bit more complicated. Indeed, using the commutator formula for horizontal derivatives (Theorem 3.5.2 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF]), we see that

| h ∇ U(x, ξ)| ≤ Cε | ξ, ν | for (x, ξ) ∈ ∂ -ΩM, ν, ξ = 0. ( 5
h ∇i (HU) j k = h ∇i (ξ p h ∇p U j k ) = ξ p h ∇i h ∇p U j k = = ξ p h ∇p h ∇i U j k + ξ p -R l qip ξ q v ∇l U j k + R j lip U l k -R l kip U j l .
This can be written as

h ∇ HU = H h ∇ U + R 1 [ v ∇ U] + R 2 [U]
(5.18) with some algebraic operators R 1 and R 2 on semibasic tensors which are determined by the curvature tensor. The operator R 2 satisfies R 2 [E] = 0. Therefore the first of estimates (5.12), which is already proved, implies the inequality

|R 2 [U]| = |R 2 [U -E]| ≤ Cε (5.19)
with some constant C depending on the curvature bound. Using commutator formulas (5.17 (5.23) The first summand of the expression in braces is equal to zero by (5.15). In virtue of (5.14), the second summand of this expression is estimated as This implies tr (π ξ f (x + tξ)π ξ ) = -f (x + tξ)ξ, ξ .

Integrating the last equality with respect to t and recalling definitions (6.1) and (6.2) of S and I, we see that S[ f ] = -I[ f ]. (6.7) Substituting (6.6)-(6.7) into (6.5), we obtain

S[f ] = -I[ f -2λE].
Thus, f is in the kernel of S if and only if

I[ f -2λE] = 0. (6.8)
Now, we apply Theorem 2.2.1 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF] which states that (6.8) is equivalent to the existence of a vector field v such that f -2λE = dv, (6.9)

where d : C ∞ (R 3 ; R 3 ) → C ∞ (R 3 ; M(3)), (dv) ij = 1 2
∂v i ∂x j + ∂v j ∂x i
is the inner derivative. Theorem 2.2.1 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF] is formulated and proved for compactly supported tensor fields. Nevertheless, the same proof works for f -2λE ∈ S(R 3 ; M(3)) and gives the vector field v belonging to the Schwartz space S(R 3 ; R 3 ). Let us express λ and f through f . Applying the trace operator to the first of equations (6.4) and taking the second one into account, we see that λ = 1 3 tr f. (6.10) From (6.4) and (6.10), fij = f ij -1 3 tr f • δ ij . (6.11) Substitute (6.10)-(6.11) into (6.9) to obtain dv = f -tr f • E. (6.12) Equation (6.12) represents the overdetermined system of six first order partial differential equations in three unknowns (v 1 , v 2 , v 3 ). The solvability condition for the system is presented by Theorem 2.2.2 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF]: equation (6.12) is solvable if and only if the righthand side of (6.12) belongs to the kernel of the Saint-Venant operator. Here, we prefer to use the version R of the Saint-Venant operator which is defined by the equation before formula (2.4.6) of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF]. So, the solvability condition for (6.12) is Rh = 0, where h = f -tr f • E.

(6.13)

The operator R is defined by the formula 4(Rh) ijkl = h ik;jl -h jk;il -h il;jk + h jl;ik .

It possesses the following symmetries:

(Rh) ijkl = -(Rh) jikl = -(Rh) ijlk = (Rh) klij .

. 6 ) 3 . 1

 631 Theorem For any n ≥ 4, there exist positive numbers δ = δ(n) and ε = ε(n) such that, for any n-dimensional CNTM (M, g) satisfying k(M, g) < δ(3.7)

  .33) while the tensor h ∇ u = ( h ∇i u jk ) has the mentioned property only in the last two indices. The summands on the right-hand side of (3.32) are orthogonal to each other, so

  .47) holds for β = 1, ε = δ = 0. By continuity and by estimates |P ξ f | ≤ |f |, |π ξ f ξ| ≤ |f |, |π ξ f * ξ| ≤ |f | for |ξ| = 1, the same inequality holds for some positive c, ε, δ, and β independent of f and satisfying 1 -β -Cε -Cδ ≥ 0. Combining (3.46) and (3.47), we obtain

. 4 )

 4 This gives us the following non-uniqueness in the linear problem. If the function λ vanishes on the boundary, λ| ∂M = 0, then u(x, ξ) = λL ξ /|ξ| 2 solves the boundary value problemHu = P ξ f, u ∂ -ΩM = 0 (4.5)with f = L ∇λ and satisfiesF [f ] = u| ∂ + ΩM = 0. (4.6)The boundary value problem (4.5) coincides with (3.2) for the unit weights p = q = E, and (4.6) means that f cannot be recovered fromF [f ].The same arguments give us some non-uniqueness in the nonlinear problem. Fix a function λ ∈ C ∞ (M) vanishing on the boundary, λ| ∂M = 0, and, for (x, ξ) ∈ ΩM, define the linear operator U(x, ξ) on T Cx M whose matrix in a positive orthonormal basis (e 1 , e 2 , e 3 = ξ) is   cos λ(x) -sin λ(x) 0 sin λ(x) cos λ(x)

  23) and β is an arbitrary positive number. The last term on the right-hand side of (4.22) can be estimated by Cε|f | 2 /|ξ| 2 as follows from (3.5) and (3.10). Estimating the first term on the right-hand side of (4.22) by (3.35), we obtain 2 Reh

1 1

 1 1 and 4.2 to obtain the statement of Theorem 5.1. Condition (3.1) is satisfied by (5.7). The weight p = U -solves the boundary value problem Hp = -p(P ξ f 1 ), p| ∂ -ΩM = E which is very similar to (5.2). Theorem 5.1 is thus reduced to the following Lemma 5.2 Let (M, g) be a CNTM and a tensor field f ∈ C ∞ (τ 1 1 M) satisfy (5.3)-(5.4). Then the solution U to the boundary value problem (5.2) satisfies the estimates |U -E| < Cε, | v ∇ U| < Cε on ΩM (5.12)

Lemma 5 . 3 0 |f (τ )|dτ exp t 0 |A

 5300 Let y = (y 1 (t), . . . , y N (t)) be the solution to the initial value problemdy dt = A(t)y + f (t), y(0) = y 0 ,where f (t) is an N-dimensional vector and A(t) is an N × N-matrix. Then |y(t)| ≤ |y 0 | + t (τ )|dτ, where |A(τ )| is the operator norm of the matrix A(τ ) defined with the help of the standard norm | • | on C N . Let us adjust this statement to our geometric setting. Lemma 5.4 Let (M, g) be a CNTM and f ∈ C ∞ (β 0 m M; ΩM), A ∈ C ∞ (β m m M; ΩM) be two semibasic tensor fields. Let y ∈ C(β 0 m M; ΩM) be the solution to the boundary value problem Hy = Ay + f, y| ∂ -ΩM = y 0 with some y 0 ∈ C(β 0 m

  t))|dt for (x, ξ) ∈ ΩM, where γ = γ x,ξ . Here the norm | • | is defined in Section 2, and the constant C depends on m and (M, g).

  be estimated together but h ∇ U is unbounded near Ω(∂M).We start with estimatingh ∇ U on ∂ -ΩM. To this end we consider equation (5.2) at a boundary point (x, ξ) ∈ ∂ -ΩM. Because of the boundary condition U| ∂ -ΩM = E, equation (5.2) gives (HU)| ∂ -ΩM = P ξ f. (5.13) Let us choose boundary normal coordinates (x 1 , . . . , x n ) in a neighborhood V of the boundary point such that the boundary is determined by the equation x n = 0, x n ≥ 0 in V , and g in = δ in . Because of the boundary condition U i j | x n =0 = δ i j , the equalities h ∇α U i j | x n =0 = 0 (0 ≤ α ≤ n -1) hold, and equation (5.13) becomes (ξ n h ∇n

∇∇

  .14) v ∇ U vanishes on ∂ -ΩM as follows from the condition U| ∂ -ΩM = E, v ∇ U| ∂ -ΩM = 0to equation (5.2), we obtain v ∇ HU = (P ξ f ) v ∇ U + ( v ∇ P ξ )f U, h ∇ HU = (P ξ f ) h ∇ U + (P ξ ∇f )U.(5.16)We have used that v ∇ f = 0 and h ∇ f = ∇f since f is independent of ξ. Operators v ∇ and H satisfy the commutator formula v from the definition H = ξ i h ∇i and from the fact that v ∇ and h ∇ commute. The commutator formula for h

∇

  ) and (5.18), we write(5.16) as P ξ )f U, G = (P ξ ∇f )U -R 2 [U]. (5.22)We first consider (5.20)-(5.21) as a linear system of ordinary differential equations in coordinates of v ∇ U and h ∇ U with free terms F and G. Applying Lemma 5.4 to the system, we obtain the estimate| h ∇ U(x, ξ)| ≤ C | v ∇ U(γ(τ -(x, ξ)), γ(τ -(x, ξ)))| + | h ∇ U(γ(τ -(x, ξ)), γ(τ -(x, ξ)))|+ + 0 τ -(x,ξ) (|F (γ(t), γ(t))|+|G(γ(t), γ(t))|)dt exp C 0 τ -(x,ξ)(|R 1 ((γ(t), γ(t))|+|P γ(t) f (γ(t))|+|E|)dt .

|h∇ 3 )

 3 U(γ(τ -(x, ξ)), γ(τ -(x, ξ)))| ≤ Cε | γ(τ -(x, ξ)), ν(γ(τ -(x, ξ)) | .Let us recall the definition of the ray transformI : S(R 3 ; M(3)) → C ∞ (T S 2 ), I[f ](x, ξ) = ∞ -∞ f (x + tξ)ξ, ξ dt = ∞ -∞ f ij (x + tξ)ξ i ξ j dt for (x, ξ) ∈ T S 2 . (6.2)See Chapter 2 of[START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF] for the theory of the ray transform on the Euclidean space. Our approach to studying the kernel of S is based on the following observation:S[f ] = -I[f ]for a symmetric matrix function f with zero trace, see equation (6.7) below.Theorem 6.1 A symmetric matrix function f = (f ij (x)) ∈ S(R 3 ; M(3)) belongs to the kernel of the operator S if and only if it satisfies the system of partial differential equations R 1 [f ] := f 11;11 + 2f 12;12 + f 22;22 + f 33;11 + f 33;22 = 0 R 2 [f ] := f 11;11 + 2f 13;13 + f 22;11 + f 22;33 + f 33;33 = 0 R 3 [f ] := f 11;22 + f 11;33 + 2f 23;23 + f 22;22 + f 33;33 = 0 where the tensor notations for partial derivatives f ij;kl = ∂ 2 f ij /∂x k ∂x l are used for brevity. The theorem gives the full system {R i [f ] | i = 1, 2, 3} of local linear functionals that can be recovered from the data S[f ].Proof. Represent the matrix f asf = f + λE, tr f = 0,(6.4)where E is the unit matrix and λ ∈ S(R 3 ) is a scalar function. Then S[f ] = S[ f] + S[λE]. (6.5) As follows from definition (6.1) of the operator S, S[λE] = 2I[λE], + tξ)dt for (x, ξ) ∈ T S 2 is the ray transform of the matrix function λE. If (ξ, η, ζ) is an orthonormal basis of R 3 , then 0 = tr f = f (x)ξ, ξ + f(x)η, η + f (x)ζ, ζ and tr (π ξ f (x)π ξ ) = f(x)η, η + f (x)ζ, ζ .

  ∂/∂x in the case of M ⊂ R n with the standard Euclidean metric and of Cartesian coordinates.

						.ir j 1 ...js +
			r				s
		+	a=1	Γ ia kp u	i 1 ...i a-1 pi a+1 ...ir j 1 ...js	-	a=1	Γ p kja u i 1 ...ir j 1 ...j a-1 pj a+1 ...js ,
	where Γ i jk are Christoffel symbols. See Sections 3.4-3.6 of [5] for properties of these
	operators. Note that	v ∇ = ∂/∂ξ and

h ∇ =

The second author was supported by RFBR grant 05-01-00611 and DFG -RFBR grant 04-01-04000. The work was started when the author visited Université de Nantes June 2006, and the work was mostly done when the author visited Universität des Saarlandes in the scope of the latter grant. The author is grateful to the whole faculty of the Applied Mathematics Department of the university and especially to his landlady, Frau Reinhardt, for her hospitality.

By Lemma 4.1.2 of [START_REF] Sharafutdinov | Integral Geometry of Tensor Fields[END_REF], the estimate

holds. Combining two last estimates, we obtain

As is seen from (5.22), (5.19), and (5.4), the integral inside the braces in (5.23) can be estimated by Cε. Finally, the integral under the exponent on (5.23) is estimated by some constant. Therefore (5.23) implies the estimate

for (x, ξ) ∈ ΩM.

(5.24)

Now, we consider (5.20) as a linear system of ordinary differential equations in coordinates of v ∇ U with the free term -h ∇ U + F . Applying Lemma 5.4 to the system and using the homogeneous initial condition (5.15), we obtain the estimate

(5.25) In virtue of (5.24), the first integral on the right-hand side of (5.25) can be estimated by Cε. Estimating then the second integral by (5.4), we obtain the second of inequalities (5.12). The lemma is proved.

Kernel of the operator S

In this section, we restrict ourselves to considering symmetric matrix functions on the whole of R 3 endowed with the standard scalar product •, • .

Let M(3) be the space of complex-valued symmetric 3 × 3-matrices. Such a matrix f ∈ M(3) is considered as the linear operator f : C 3 → C 3 . By S 2 we denote the unit sphere in R 3 and by T S 2 = {(x, ξ) | ξ ∈ S 2 , x ∈ R 3 , x, ξ = 0}, the tangent bundle of the sphere. Given ξ ∈ S 2 , let ξ ⊥ = {η ∈ C 3 | η, ξ = 0} be the complex two-dimensional space of vectors orthogonal to ξ and π ξ : C 3 → C 3 , the orthogonal projection onto ξ ⊥ .

Let S(R 3 ; M(3)) be the Schwartz space of M(3)-valued functions on R 3 . The linear operator S :

is defined by

compare with (1.10). We are going to answer the question: to which extent is a symmetric matrix function f ∈ S(R 3 ; M(3)) determined by the data S[f ]? Since S[f ] depends linearly on f , the question is equivalent to studying the kernel of the operator S.

Because of the symmetries, the tensor Rh has six linearly independent components, and equation ( 6.13) is equivalent to the system

)), system (6.14) contains only three independent equations. More precisely: each of the last three equations of (6.14) can be obtained from the first three equations by taking linear combinations, differentiation, and integration. To prove this, let us rewrite system (6.14) in terms of the Fourier transform g(ξ) = f . Applying the Fourier transform to each equation of (6.14), we arrive to the system

One can easily see the following three relations between equations of system (6.15):

Therefore three last equations of (6.15) follow from three first equations at least if g(ξ) depends continuously on ξ.

Deleting three last equations from system (6.15), we obtain the equivalent system

The same is true for system (6.14): deleting three last equations from (6.14), we will obtain the equivalent system (6.3). The theorem is proved.

System (6.16) enables us to answer the question: which integral moments of f can be determined from S[f ]? Indeed, the Tailor series of the function g(ξ) = f is

where µ

are the integral moments of order m. Assuming f to be in the kernel of S, let us insert series (6.17) into system (6.16). Since the coefficients of (6.16) are homogeneous functions of ξ, the system does not mix moments of different orders. This means that we can take g(ξ) in the form

if we are looking for moments of order m.

Let us start with considering zero order moments. We substitute the expressions (6.16). Equating coefficients at the same degrees of ξ at the resulting equations, we easily find that µ (0) ij = 0 for every (i, j). This means that the integral R 3 f (x)dx can be determined from the data S[f ].

Next, we consider first order moments. We substitute the expressions

ij,3 ξ 3 into system (6.16). Equating coefficients at the same degrees of ξ at the resulting equations, we arrive to the system

22,1 = 0, µ

11,1 + µ

33,1 = 0, µ

11,2 + µ

22,2 = 0, µ

22,2 + µ

33,2 = 0, µ

11,3 + µ

(1)

33,3 = 0, µ

22,3 + µ

11,2 + 2µ

(1)

12,1 + µ

33,2 = 0, µ

11,3 + 2µ

(1)

13,1 + µ

22,3 = 0, µ

11,3 + µ

12,2 + µ

22,1 + µ

33,1 = 0, 2µ

13,3 + µ

22,1 + µ

33,1 = 0, µ

11,2 + 2µ

(1)

33,2 = 0. The general solution to the system looks as follows:

13,1 = a 3 , µ

12,2 = a 1 , µ

13,2 = 0, µ

22,2 = a 2 , µ

23,2 = a 3 , µ 33,3 = a 3 , where (a 1 , a 2 , a 3 ) are arbitrary constants. Eliminating the constants, we obtain the following independent system of 15 linear combinations of first order moments of f which can be recovered from the data S[f ]: (µ ii,i )[f ], (6.18) where δ ij is the Kronecker tensor. System (6.18) is considered for such (i, j, k) that at least two of these indices are different. A similar consideration is possible for integral moments µ


Benoît Pier
email: benoit.pier@ec-lyon.fr

S Derebail Muralidhar

J F Scott

R Govindarajan

Flow around a rotating, semi-infinite cylinder in an axial stream

Keywords: fluid mechanics, applied mathematics rotating body, wall jet, boundary layer, asymptotic analysis

This paper concerns steady, high-Reynolds-number flow around a semi-infinite, rotating cylinder placed in an axial stream and uses boundary-layer type of equations which apply even when the boundary-layer thickness is comparable to the cylinder radius, as indeed it is at large enough downstream distances. At large rotation rates, it is found that a wall jet appears over a certain range of downstream locations. This jet strengthens with increasing rotation, but first strengthens then weakens as downstream distance increases, eventually disappearing, so the flow recovers a profile qualitatively similar to a classical boundary layer. The asymptotic solution at large streamwise distances is obtained as an expansion in inverse powers of the logarithm of the distance. It is found that the asymptotic radial and axial velocity components are the same as for a nonrotating cylinder, to all orders in this expansion.

Introduction

When a semi-infinite rotating cylindrical body is placed in a high-Reynolds-number axial flow (see figure 1), an axisymmetric boundary layer develops along the cylinder. Initially thin, this layer becomes of thickness comparable with, then larger than the cylinder radius at sufficiently large axial distances. Our original motivation for studying this flow was to undertake a stability analysis. However, it soon became clear that there are very few existing studies of the underlying flow, c The Author(s) Published by the Royal Society. All rights reserved.

In §2, we define the two non-dimensional control parameters of the problem, Re and S, the Reynolds number and non-dimensional rotation rate. The boundary-layer equations, valid for large Re, and allowing for boundary-layer thickness to be comparable with the cylinder radius are given. These equations generalise the Prandtl equations and apply for arbitrary (not necessarily small) ratios of boundary-layer thickness to cylinder radius. Suitable rescaling of the variables renders the problem independent of Re, leaving only S as control parameter. Section 3 describes the numerical scheme and its verification, while §4 gives results and discussion, in particular focussing on the wall jet. Finally, §5 gives asymptotic analyses of the limits of large Z = z/Re and large S. The boundary layer on the nose is discussed in the appendices. It is found that the precise shape of the nose is unimportant: the input to the boundary-layer equations of §2 being the Blasius flat-plate flow (generalized to include the azimuthal component due to rotation), independent of the nose shape.

Boundary-layer equations

A semi-infinite cylinder of radius a rotates about its axis with angular velocity Ω and is placed in an axial stream of incompressible fluid of velocity U∞ and viscosity ν (see figure 1). Assuming large Reynolds number, an initially thin boundary layer develops along the cylinder. If the cylinder were sharply truncated at the nose, flow separation would occur as is usually the case at a salient edge [START_REF] Batchelor | An Introduction to Fluid Dynamics[END_REF]. To avoid this scenario we assume there is a smooth nose, as shown in figure 1. Another way to avoid flow separation would be to consider a thin hollow cylinder. The boundarylayer equations formulated in this section hold good for both these cases. A detailed analysis of the nose region is given in appendix B. The natural length and velocity scales are a and U∞. These scales are used to nondimensionalize the axisymmetric, steady Navier-Stokes equations in cylindrical coordinates, z, r, θ. There are two non-dimensional parameters, namely the Reynolds number

Re = U∞a ν , (2.1)
and the rotation rate

S = Ωa U∞ . (2.2)
Assuming a large Reynolds number, the length scale for axial variation of the flow is much longer than that for radial variation. This separation of scales leads to the boundary-layer approximation. Thus,

uz ∂uz ∂z + ur ∂uz ∂r = - ∂p ∂z + 1 Re ∂ 2 uz ∂r 2 + 1 r ∂uz ∂r , (2.3)
u 2 θ r = ∂p ∂r ,
(2.4)

uz ∂u θ ∂z + ur ∂u θ ∂r + u θ ur r = 1 Re ∂ 2 u θ ∂r 2 + 1 r ∂u θ ∂r - u θ r 2 , (2.5
)

∂uz ∂z + 1 r ∂rur ∂r = 0 (2.6)
are obtained by dropping terms of higher order from the Navier-Stokes equations in the usual manner. Note that we have not assumed the boundary layer to be thin compared with the radius. The boundary conditions are

uz = 0, ur = 0, u θ = S z > 0, r = 1, (2.7
)

uz → 1, u θ → 0, p → 0 z > 0, r → ∞. (2.8)
The above equations contain the azimuthal component, u θ , of the velocity. This is due to rotation of the cylinder, which induces the centrifugal term on the left-hand side of equation (2.4), leading to a significant radial pressure gradient. Such an effect is not present in classical boundary-layer theory, which predicts near constancy of the pressure across the layer. Compared to the Prandtl equations of a classical boundary layer, equations (2.3)-(2.6) allow for the additional effects of both rotation and curvature. Near the nose the boundary layer is thin compared with the cylinder radius and curvature effects are negligible. But at large z, boundary-thickening eventually makes the thickness comparable to, then larger compared with the radius, and the full set of equations is required.

The above problem needs to be completed by inlet conditions. At distances from the nose of O(a), the boundary layer is thin compared to the radius and is described by the axisymmetric Prandtl equations given in appendix B. On the cylinder (after leaving the nose), these equations become the flat-plate Prandtl equations, and as z increases, we expect the flow to forget the precise initial conditions and to approach the Blasius solution, independent of the nose shape (here, we implicitly suppose the nose length to be of the same order as its diameter). There are, in fact, two asymptotic regions, z = O(1), where the equations of appendix B apply, and z = O(Re), where the boundary-layer thickness is comparable to the cylinder radius and equations (2.3)-(2.8) hold. Matching between these regions requires that the inlet condition be the Blasius solution (generalized to include the azimuthal component due to rotation). The same is true for the case of the hollow cylinder. In either case, equations (2.3)-(2.8) are supplemented by Blasius initial conditions as z → 0.

Introducing the scaled variables

R = r, Z = z Re , (2.9
)

Uz = uz, Ur = RurRe, U θ = Ru θ , P = R 2 p, (2.10
)

equations (2.3)-(2.8) become Uz ∂Uz ∂Z + Ur R ∂Uz ∂R = - 1 R 2 ∂P ∂Z + ∂ 2 Uz ∂R 2 + 1 R ∂Uz ∂R , (2.11)
U 2 θ = R ∂P ∂R -2P, (2.12)
Uz ∂U θ ∂Z + Ur R ∂U θ ∂R = ∂ 2 U θ ∂R 2 - 1 R ∂U θ ∂R , (2.13)
∂Uz ∂Z + 1 R ∂Ur ∂R = 0, (2.14)
Uz = Ur = 0, U θ = S R = 1, (2.15) Uz → 1, U θ = 0, P = 0 R → ∞. (2.16)
It is apparent that, using these scalings, Re has disappeared from the problem, leaving S as the only non-dimensional parameter. This result indicates, among other things, that the natural scaling of the axial coordinate is z = O(Re). Thus, as noted earlier, the distance needed for the boundary-layer thickness to become comparable with the radius is Re times the radius. The factors of R appearing in equation (2.10) have been introduced to improve numerical convergence.

Numerical scheme and validation

The boundary-layer thickness goes to zero like Z 1/2 and Ur → ∞ like Z -1/2 as Z = 0 is approached. To maintain numerical accuracy in the presence of such singular behaviour, we introduce the variables

ζ = (2Z) 1 2 , σ = R -1 ζ , Vr = ζUr, Vz = Uz, V θ = U θ . ( 3
ζVz ∂Vz ∂ζ + Vr -ζ R -σVz ∂Vz ∂σ = 1 R 2 (σ ∂P ∂σ -ζ ∂P ∂ζ ) + ∂ 2 Vz ∂σ 2 , (3.2) ζ R V 2 θ + 2P = ∂P ∂σ , (3.3)
ζVz ∂V θ ∂ζ + Vr + ζ R -σVz ∂V θ ∂σ = ∂ 2 V θ ∂σ 2 , (3.4)
ζ ∂Vz ∂ζ -σ ∂Vz ∂σ + 1 R ∂Vr ∂σ = 0, (3.5)
with the boundary conditions

Vz = Vr = 0, V θ = S σ = 0, (3.6
)

Vz → 1, V θ = 0, P = 0 σ → ∞. (3.7)
These equations govern the axial evolution of the flow. The inlet condition (Blasius solution) is obtained by setting ζ = 0 and solving the resulting equations.

The radial coordinate σ is discretized using Chebyshev collocation points:

xn = cos nπ N -1 0 ≤ n < N, (3.8
)

σn = σ(1 + xn) 1 -xn , x ∈ [-1, 1] → σ ∈ [0, ∞]. (3.9)
The parameter σ controls the distribution of points such that half of them lie between 0 ≤ σ ≤ σ.

The velocities Vz and V θ are represented by their values at all collocation points. However, since there is no boundary condition for the pressure at the surface, it is represented at all points except σ = 0. Similarly, there is no boundary condition for Vr at σ = ∞ and so it is represented at all points apart from σ = ∞. The Chebyshev derivative matrices for P and Vr are correspondingly modified (e.g. appendix A in [START_REF] Scott | Moffatt-type flows in a trihedral cone[END_REF]).

The coordinate ζ is discretised using small, equally-spaced steps, ζ i = i∆, and the variables Vz, Vr, V θ , P are represented by their values at ζ i . Equations (3.2), (3.4) and (3.5) are evaluated at mid-step, ζ i+ 1 2 , using an implicit scheme that employs centered finite differencing to represent the ζ-derivatives. Equation (3.3) is evaluated at the step position ζ i , rather than at the midstep. At each step, the equations are solved using Newton-Raphson iteration, thus allowing forward marching. The inlet solution is obtained from equations (3.2)-(3.7) using ζ = 0. Following discretization in σ using the collocation points, the result is again obtained by Newton-Raphson iteration.

The code was first tested by changing the numerical parameters N, σ, ∆, and observing the dependence of the solution on these parameters. Based on the convergence results we decided to use N = 128, σ = 5, ∆ = 0.001 in our computations. These values gave convergence to better than seven decimal places. The code was also tested using the volume-flux and momentum balance equations. The results respect these equations to seven decimal places. Although use of the Blasius solution at the inlet has earlier been justified by an asymptotic argument, it is interesting to see the effect of a change in inlet profile on the solution. Thus we modified the inlet profile to be U * z = Uz + Aσexp(-0.5σ), where Uz is the Blasius profile. Taking A = 2, it was found that the change in the Uz at Z = 2 was in the sixth decimal place. This illustrates the fact that the flow forgets the initial condition as Z increases and becomes insensitive to the precise inlet profile used.

We also validated the code by comparing our results with the existing literature. Tutty et al. [START_REF] Tutty | Boundary layer flow on a long thin cylinder[END_REF] studied the case without rotation. The axial (x t ) and radial (σ t ) coordinates used by Tutty et al. [START_REF] Tutty | Boundary layer flow on a long thin cylinder[END_REF] are related to those used here via x t /Re = Z and σ t = √ 2σ. Figure 2 shows good agreement with our results for Re = 10 4 , x t = 0.01 and x t = 10 5 . Herrada et al. [START_REF] Herrada | Stability of the boundary layer flow on a long thin rotating cylinder[END_REF] considered the problem of the rotating cylinder. They do not give velocity profiles, but rather the skin friction on the cylinder:

τ = ∂Uz(Z) ∂R R=1 = 1 ζ ∂Vz(ζ) ∂σ σ=0 . ( 3

Results

Flow profiles were obtained for different values of S and Z. Figure 3 shows velocity profiles for S = 0.1 and different values of Z. For Z = 0, Vz(σ) corresponds to the Blasius flat-plate solution.

As Z increases, Vz(σ) deviates from the Blasius profile due to cylinder curvature and rotation. It should be borne in mind that the boundary-layer thickness increases with Z, although this is not apparent in the figure because the scaled radial coordinate σ = (R -1)/ζ has been used. Note that the azimuthal velocity at Z = 0 is V θ (σ) = S 1 -Vz(σ) . As Z increases, small departures from this profile arise. Figure 4 and figure 5 show results for Z = 0.5 and different values of S.

When S 1, Vz(σ) is a modified Blasius profile. However, for S 4, Vz(σ) is no longer monotonic having a maximum at finite σ. At large S, the maximum is large and the profile is better described as an axial wall jet, rather than a boundary layer. At first sight, it is perhaps surprising that increasing the rotation rate leads to a stronger and stronger axial flow. Increasing S causes V θ to increase (see figure 4). This in turn produces an increasing radial pressure gradient due to the centrifugal force. Since the pressure is constant outside the boundary layer, the pressure within the layer drops (see figure 5) with Z. The development of the flow means that the axial presure gradient becomes larger and larger, thus driving a strong axial wall jet. Whereas for lower values of S, boundary-layer thickening leads to positive Vr, at large S entrainment by the wall jet gives negative Vr outside the layer (see figure 4).

The existence of a wall jet at large S is apparent in the axial velocity profiles given by Petrov [START_REF] Petrov | Boundary layer on a rotating cylinder in axial flow[END_REF]. However, that article only gives such profiles for the case of small Z (thin boundary layer) and no explanation is provided. Petrov [START_REF] Petrov | Boundary layer on a rotating cylinder in axial flow[END_REF] also gives the maximum axial velocity U max z = max R (Uz(R)) as a function of Z. Figure 6 shows a comparison with our results. A small difference is apparent, the origin of which is unclear.

Figure 7 shows contours of constant U max z in the (S, ζ)-plane as well as the boundary (solid line) separating the region in which U max z = 1 from that in which U max z > 1 (which we interpret as indicating a wall jet). It will be seen that there is a threshold, S = 4.15, below which U max z = 1. Above this value, the wall jet exists for some range of axial position. Note that, whatever the strength of rotation, the wall jet eventually disappears sufficiently far downstream.

The thickness of the boundary layer/wall jet can be measured using

δ = 1 U max z ∞ 1 |1 -Uz|dR. (4.1)
The absolute value is taken to make the integral always positive and the division by U max z allows for the strong wall jets which arise at large S. Figure 8 shows δ(Z) for different values of S. The layer thickness is seen to increase with Z in a roughly parabolic manner (recall that the thickness behaves as Z the usual manner. Increasing S causes the layer to become thinner. At large S, the wall jet is of increasing strength. Viscous diffusion competes with axial convection, the latter being of growing importance, hence the decrease of δ with increasing S. the streamwise direction as Z increases.

Asymptotic analysis (a) Large-Z asymptotics

U 2 θ + 2P = η ∂P ∂η , (5.3)
Uz ∂U θ ∂χ -η ∂U θ ∂η + Ur η ∂U θ ∂η = ∂ 2 U θ ∂η 2 - 1 η ∂U θ ∂η , (5.4
Uz = Ur = 0, U θ = S η = e -χ , (5.6
)

Uz → 1, U θ → 0, P → 0 η → ∞. (5.7)
It is shown in appendix A that Uz, Ur, U θ and P have asymptotic expansions in powers of χ -1 . The factor of e -2χ in equation (5.2) is exponentially small and is hence negligible at all algebraic orders. Without the corresponding term in equation (5.2), Uz and Ur decouple from U θ and P , though the latter depends on the former. Thus, we expect such decoupling to hold at all orders. This is indeed what is found in appendix A, where the governing equations for the coefficients of the expansions in powers of χ -1 are obtained for all orders. Given decoupling, rotation does not enter into the asymptotics of Uz and Ur, which are consequently the same as for the non-rotating case.

Glauert & Lighthill [START_REF] Glauert | The axisymmetric boundary layer on a long thin cylinder[END_REF] studied the case without rotation and obtained the expansions of Uz and Ur. Appendix A extends the analysis to include rotation and gives detailed results up to order 5. At first order, the asymptotic solution can be obtained analytically and is given by

Uz ∼ 1 -χ -1 ∞ η e -ξ 2 /2 ξ dξ, (5.8
)

Ur ∼ χ -1 (1 -e -η 2 /2 ), (5.9)
U θ ∼ Se -η 2 /2 , (5.10) .11) In figure 9, the results for Uz show convergence to the asymptotic form (5.8), while those for u θ = U θ /R converge to u θ ∼ S/R, which is the flow due to a rotating cylinder, infinite in both axial directions (rather than semi-infinite) and without axial flow.

P ∼ -S 2 η 2 ∞ η e -ξ 2 ξ 3 dξ. ( 5

(b) Large-S asymptotics

As we saw in the previous section, the numerical results show the existence of a wall jet at large S.

In this limit, appropriate scaled variables are

Z * = Z S , R * = R, (5.12
) order (S 0 ), we find

U * z = Uz S , U * r = Ur, U * θ = U θ S , P * = P S 2 . ( 5
U * z ∂U * z ∂Z * + U * r R * ∂U * z ∂R * = - 1 R * 2 ∂P * ∂Z * + ∂ 2 U * z ∂R * 2 + 1 R * ∂U * z ∂R * , (5.14)
U * 2 θ = R * ∂P * ∂R * -2P * , (5.15)
U * z ∂U * θ ∂Z * + U * r R * ∂U * θ ∂R * = ∂ 2 U * θ ∂R * 2 - 1 R * ∂U * θ ∂R * , (5.16
)

∂U * z ∂Z * + 1 R * ∂U * r ∂R * = 0, (5.17)
with the following inlet and boundary conditions

U * z = 0, U * θ = 0 Z * = 0, (5.18)
U * z = 0, U * r = 0, U * θ = 1 R * = 1, (5.19)
U * z = 0, U * θ = 0, P * = 0 R * → ∞.
(5.20)

Figure 10 shows the solution of the above problem (solid line) compared with the numerical results discussed before for Z * = 0.1 and different values of S. It is apparent that the asymptotics are indeed approached as S → ∞. Figure 11 shows the leading-order asymptotic solution for different values of Z * . We see that the large-Z limit (u θ ∼ SR -1 ) is approached by u θ as Z * → ∞.

Note that the limit S → ∞ can be reached in two ways: either by increasing the rotation rate, or by decreasing the velocity U∞ to zero. Note also that uz/S and u θ /S are the velocity components non-dimensionalized by Ωa, rather than U∞, and that Z * = z/Re Ω , where Re Ω = Ωa 2 /ν is the Reynolds number based on the rotational velocity Ωa. Thus figure 11 can be interpreted as showing the flow due to a rotating, semi-infinite cylinder in a still fluid (U∞ = 0). It can be shown that the separation of radial and axial length scales, which underlies the boundary-layer type approximation we have used, is valid if either of the Reynolds numbers, Re or Re Ω , is large.

Conclusion

In this paper we have presented a study of the flow around a rotating cylinder in an axial stream. We have assumed a smooth nose to avoid flow separation. The two non-dimensional control parameters of the problem are: Reynolds number (Re) and rotation rate (S). The flow equations are formulated using a boundary-layer type approximation, appropriate at large Reynolds numbers and in which the flow is assumed to evolve slowly in the streamwise direction in comparison to the radial direction. The resulting equations are not limited to the case in which the boundary layer is thin compared with the cylinder radius. By using appropriate scalings we remove Re from the problem.

The results show that the boundary-layer thickness increases with axial distance, becoming comparable with the cylinder radius a at distances of O(Re a). Prior to this, the layer is thin compared to the radius and the flow is close to the Blasius profile of a flat plate. However, it differs from the Blasius solution due to effects of curvature and rotation at larger downstream distances. As S increases, the centrifugal force creates an increasing radial pressure gradient, which combined with axial development, implies an increasing axial gradient of pressure. Above S = 4.15, the maximum velocity exceeds the free-stream velocity for a range of Z and we say that a wall jet exists. This jet becomes stronger and stronger as S → ∞.

In the limit of large Z, we find that the axial and radial components of velocity decouple from the azimuthal velocity component and pressure. All these quantities are found to have asymptotic expansions in inverse powers of ln(Z), a result already obtained for the non-rotating case by Glauert & Lighthill [START_REF] Glauert | The axisymmetric boundary layer on a long thin cylinder[END_REF], and here extended to include rotation. The leading-order term in the u θ expansion is u θ ∼ SR -1 , which is the flow expected for a rotating cylinder, infinite in both directions. Because U θ becomes independent of Z, the same is true of the pressure field resulting from the centrifugal force, hence the absence of an axial pressure gradient to drive the axial/radial flow. This is the reason for the decoupling.

When S is large, we introduce appropriate scalings for Z, Uz, U θ and P . The asymptotic expansions of the scaled velocity and pressure proceed as inverse powers of S, beginning with S 0 . The leading-order term describes an axial wall jet due to a rotating cylinder in a fluid at rest.

(n) z , U (n) r , U (n) θ
and P (n) . A 1 = 1 gets the process started and leads to the leading-order outer solution, (5.8)- (5.11), in agreement with [START_REF] Glauert | The axisymmetric boundary layer on a long thin cylinder[END_REF]. Equation (A.29) gives An at the next order, allowing solution at successively higher orders. It can be shown that

U (n) z = Anlnη + A n+1 + O(η 2 ln pn η), U (n) θ = 1 -qn + O(η 2 ), U (n) r = O(η 2 ln qn η), (A.30)
as η → 0, where qn = p n+1 -1, p 1 = 0, p 2 = 1 and pn = 2 for n ≥ 3. The terms in (A.30) indicated by the O() notation are exponentially small in the inner region, while the remaining ones reproduce the inner solution. Thus, the outer expansions in fact apply in the inner region.

The above procedure has been implemented numerically and results up to n = 5 are presented in figure 12. Figure 13 shows the comparison of numerical solution of Uz at Z = 5000 and S = 1 with the asymptotic solution obtained by truncating at different orders n. Although this result shows good convergence, and therefore further confirms both numerical and analytical results, it should be borne in mind that the expansions (A.1)-(A.4) are, in fact, asymptotic as Z → ∞, rather than necessarily convergent at any finite Z.

B. Flow over the nose

Large Reynolds number implies a thin boundary layer over the nose. Schlichting [START_REF] Schlichting | Laminar flow about a rotating body of revolution in an axial airstream[END_REF] gives the axisymmetric boundary-layer equations in terms of curvilinear coordinates, x, y, θ, where x is distance along the surface, and y is distance normal to the surface. Here we use the non-dimensional coordinates, velocity components and pressure: ỹ = Re

Figure 1 .

 1 Figure 1. Schematic diagram of the problem.

. 10 )Figure 2

 102 Figure2shows τ as a function of Z for S = 1 and good agreement is apparent.

Figure 2 .

 2 Figure 2. Comparison of Uz for S = 0 at two values of Z with [5] and comparison of skin friction τ for S = 1 with [9].

Figure 3 .

 3 Figure 3. Velocity profiles Vz(σ), V θ (σ) and Vr(σ) for different axial positions Z at S = 0.1.

1 2Figure 4 .

 14 Figure 4. Velocity profiles Vz(σ), V θ (σ) and Vr(σ) at Z = 0.5 for different values of rotation rate S.

Figure 5 .Figure 6 .

 56 Figure 5. Pressure profile P (σ) at Z = 0.5 for different values of S.

  ) rspa.royalsocietypublishing.org Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 7 .

 7 Figure 7. Grayscale plot of U max z alongside a contour plot of U max z . The figure also shows the boundary (solid line)

Figure 8 .

 8 Figure 8. Boundary-layer/wall-jet thickness δ versus Z for different values of S.

Figure 9 .

 9 Figure 9. Comparison of χ(Uz -1) and u θ = U θ /R for different values of Z at S = 0.1 with the Z → ∞ solution.

  .13) The scaling of Z reflects the increasing distance required for flow development as the rotation rate increases. The scaling of Uz and U θ indicates the strengthening flow velocity as S increases. The large-S asymptotic expansions of U * z , U * r , U * θ and P * proceed as powers of S -1 . At leading rspa.royalsocietypublishing.org Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 10 .Figure 11 .

 1011 Figure 10. Comparison of U * z and U * θ for different values of S and Z * = 0.1 with the S → ∞ asymptotic solution.

Assuming

  rspa.royalsocietypublishing.org Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . as η → 0. It follows from (A.25) that η dU An is known, equations (A.5)-(A.11) and the boundary conditions (A.12) and (A.26)-(A.28) can be solved for U

1 2Figure 12 .

 112 Figure 12. First five coefficients of the large-Z asymptotic expansions of Uz, U θ and Ur.

  .1)Here, the boundary-layer thickness is prevented from going to zero in the radial coordinate σ by dividing R -1 by ζ. Vr is kept finite by use of the factor ζ and ζ is used in place of Z to make the solution a smooth function of the axial coordinate. Using these variables in the boundary-layer

	equations (2.11)-(2.16) results in

Funding statement. The authors are grateful to the French Ministry of Education and Research and the Indo-French Centre for the Promotion of Advanced Research-CEFIPRA for financial support.

A. Large Z asymptotic expansions

The flow variables are expressed as asymptotic expansions in inverse powers of χ:

Introducing these expansions into equations (5.2)-(5.5) gives

where

The boundary conditions at η → ∞ are

Application of the boundary conditions (2.15) at the cylinder surface requires the introduction of an inner region, R = O(1), represented by the expansions

The boundary conditions (2.15) imply Bn = 0, C 1 + D 1 = 1, and Cn + Dn = 0 for n > 1.

Recalling that R = ζη = e χ η, the inner expansions give