N Reuge'

F M Schmidt2

Y Le Maoulp

M Rachik3

F Abbe4

Elastomer Biaxial Characterization Using Bubble Inflation Technique. I: Expe rim enta I Investigations

fiance

The biaxial rheological behavior of materials such as elastomers or polymers can be obtained using a bubble-inflation-technique. A circular membrane is clamped at the rim and expanded under gas pressure. The inflation of the circular membrane is recorded using a CCD video camera and the blowing pressure by a pressure sensor, Then, fi-om elongation and curvature radius measurements at the pole of the bubble, one can deduce equibiaxial stress-strain data. This study describes the optimization of a bubble-inflation rheometer. The most sensitive point of the technique is the estimation of the elongation at the bubble pole, deduced from video camera measurements. A direct measurement of the bubble thickness was performed using a magnetic probe in order to validate rheometer results. Such a validation has evidently never been carried out before. Results of quasi-static equibiaxial characterization of elastomers are presented and analyzed.

1) INTRODUCTION

here is a growing interest in determining biaxial T properties of materials such as elastomers used for mechanical applications or thermoplastics for the simulation of plastic-processing operations such as blow molding, stretch blow molding or thermoforming [START_REF] Delorenzi | computers & Stmctures[END_REF]. In these specific processes, the polymer is heated above the glass-transition temperature but far from the melting point. Consequently, the biaxial rheometer should be suitable for a material close to a solidlike state. In addition, the characterization of elastomers is a difficult task because of their high *Corresponding author. 522 deformation ability and their complex behaviors (nonlinear elasticity, viscoelasticity, crystallhation, instability, and so on).

Different types of biaxial rheometers are used toward this end. The classic technique of biaxid plane extension of cruciform samples is not really suitable for such materials. The Meissner rheometer (2), made of rotary clamps, seems to be more suitable (but complex and expensive, especially in multiaxial configuration). Another way to characterize such materials is the squeeze-flow-technique, where a sheet of sample is compressed. A special lubricant is required between the sample surfaces and the compressing plates in order to avoid any shear flow within the sample. But this rheometer is more convenient for liquid-like materials.

POLYMER ENGINEERING AND SCIENCE, MARCH 2001, Vol. 41, No. 3 Ultimately, bubble inflation seems to be the best way to carry out an equibiaxial elongation of elastomers with a relatively simple device.

The inflation of circular membranes has been of interest to scientists and engineers for many years. Bubble inflation with compressed air was used for the first time by Treolar (3) in 1944, to study the behavior of natural rubber and its mechanisms of bursting. In 1951, Rvlin et al. (4) used this technique to identify Neo-Hookean and Mooney-Rvlin rheological parameters, and showed that the equibiaxial state of deformation at the pole bubble is equivalent to compression according to the normal axis of the membrane. In 1971, Alexander (5) studied and modeled atmospheric neoprene balloon instabilities using the bubble inflation technique. In 1973, Paisant (6) studied the deformation of natural rubber and SBR membranes to model their biaxial behaviors.

Polymer materials have been also investigated with this technique. In 1975, Schmidt and Carley (7) developed numerical models using non-linear constitutive equations issued from the rubber-like materials theory for the membrane-inflation problem and compared their results with experimental values obtained with high-impact polystyrene and cellulose acetate butyrate above their respective glass-transition temperatures. In 1980, Denson and Hylton (8) used oil to inflate ABS membranes near the forming temperature (138°C). Oil being incompressible, they could cany out retardation experiments, which were impossible with air. In 1992, Thomas and Tintillier (9) used a video camera system (CCD camera) and preprinted membranes to analyze equibiaxial deformations of PET membranes heated above the glass transition temperature. In 1994, Kappel(l0) used the same kind of device to study an elastomer. They reached promising results, but the equibiaxial rheometer is not yet really operational.

Q. 1. Inflated bubble in biaxial ewtenswn 9

In 1996, Michaeli et at. (1 1) developed a specific bubble inflation rheometer, in which again, a membrane is inflated using silicon oil with some improvements to prescribe a constant strain rate during the entire process.

In 1997, Verron [START_REF] Verron | Contribution experimentale et numerique au proctdes de moulage par souMage et thermoformage[END_REF] developed a numerical procedure using bubble height to compute the rheological parameters for a given constitutive equation during membrane inflation. This numerical approach does not require a complex video camera device to mark the membrane, but the method requires the knowledge of the constitutive equation.

Our goal is to develop and optimize a rheometer based on a bubble-inflation technique similar to those of Thomas, Tintillier and Kappel. For this, we have proceeded to comparisons with the measured thickness of the inflated membrane using the inflation device and a magnetic probe thickness measurement. The inflation rheometer was then applied to determine the biaxial rheological properties of elastomers at ambient temperature (no heating stage).

2) BUBBLE I"LATI0N ASSulvIpTIONS

2.1) Hoop Stresa Calculation

A uniform circular membrane is clamped at the rim and inflated using compressed air applied to one side in order to produce a bubble as sketched in Rg. 1. The inflation of the circular membrane results in a truly biaxial stretching at the pole of the bubble and a planar extension near the rim. The assumption of spherical symmetry allows us to equate a,, to uLPq at the pole of the bubble. So, the Cauchy stress tensor a_expressed in spherical coordinates reduces to:

The ratio between the thiclmess of the inflated membrane e and the curvature radius R, being small enough [START_REF] Schmidt | [END_REF]14), the thin shell assumption allows us to neglect the radial stress u , in front of the hoop stress ow. In addition, we equate uee to the thicknessaverage hoop stress, which leads to: assuming that the thickness is uniform at the pole of the bubble. AP is the differential inflation pressure. uo is an additional stress due to the initial compression of the membrane.

2.2) Thioknerrs chhdation

The differential inflation pressure may be easily recorded using a pressure sensor and the curvature radius measured using video camera images analysis in order to calculate the hoop stress (Eq 2). For the thickness calculation, we need to locate a marked area around the pole of the bubble where the deformation of the membrane is equibiaxial and the thickness uniform (See Section 3).

Using the volume balance at a local level, the thickness may be computed using the area change:

(3) where e, is the flat circular membrane thickness. The hoop stretch A , , is the ratio between the current length of the profile contour of the polar zone I and the initial length of the polar zone 6: e0 e=-

GI (4)

1 he, = < Then, the thickness-average hoop stress leads to:

(5)

3) BUBBLE INFLATION TECHNIQUE

A schematic view of the bubble inflation rheometer is presented in Fig. 2. A circular membrane is mounted between two metal disks containing a circular aperture (40 mm diameter) and clamped on a support. The membrane has a 60-mm diameter and a 2-mm thickness. The air pressure is regulated with a very accurate device and recorded using a pressure sensor located close to the backside of the membrane. The inflation of the membrane is recorded using a CCD SONY video camera (752 X 582 pixels). A PC computer is used to send a s@ to the pressure regulator and to acquire images of video camera by a MATROX image grabber (24 frames/s). Elastomeric membranes were successfully inflated beyond the hemispherical stage to large elongation up to 4 (&.

3). It appears immediately that the hemispherical deformation occurs only at low elongation (5 2.5)

A very sensitive point is that elongation and curvature radius must not be measured on the entire bubble contour because deformation of the membrane is rigorously equibiaxial only at the pole of the bubble. Consequently, we had to mark the bubble in order to locate a polar area where equibiaxiality can be considered (this assumption will be verified in Section 4). A viscous silicon paste was used to perform a rectangular white mark with an initial width of 10 mm. So, the polar area could be extracted, the length of its contour measured, and its curvature radius computed using a least squares fit (Q. 4). These operations were carried out by our image-processing software (15). Finally, the hoop stretch and stress could be deduced from videometric and pressure measurements.

4) BIAXIAL REEOMETER ASSESSMENT

An experimental setup was carried out in order to validate the polar area assumption. For this, we have to compare the thickness measured at the pole of the bubble using the bubble-inflation rheometer (video camera measurement) and the thickness measured directly by another device, using a magnetic probe. Such a validation seems never to have been carried out before.

The standard technique used for the thickness measurement requires a magnet ball and a magnetic probe [START_REF] Lejeune | Contrijle dtpaisseur[END_REF]. The magnet ball is placed upon the membrane and the clamping system is reversed (Fig. 5). The magnetic field measured by the probe depends only on the distance between the probe and the ball, which represents the thickness of the membrane. For a 2-mmthick sample, the uncertainty of this device is 530 pm. Then, the membrane is inflated for increasing pressure values. For each pressure value, the magnetic probe fixed on a rack is in perfect contact with the pole of the reversed bubble. We proceed to five measurements for each pressure value. The pole bubble height H is measured using the video camera image analysis.

4.1) Spherical Bubble-Inflation Amsumption

If we assume that the bubble remains spherical during inflation (the curvature radius and the thickness are assumed to be constant), the volume balance on the entire bubble (aesumption 1) permits calculation of the relative thickness -: In Rg. 6, the relative thicknessmeasured using the bubble-inflation rheometer and using assumptions 1 and 2 are plotted versus bubble height. In both cases, the error between the measured thickness and the calculated thickness is important. The maximum error be is more than lOO?!. Thus, the spherical bubble ase sumption together with the volume balance assumption will not yield useful values for this application. Note that for the assumption 2, the discrepancy is lower than for the assumption 1.

4.2) Polar Area Assumption

In Fig. 7, the relative thickness is plotted versus bubble height for the magnetic probe thickness and the rheometer thickness. Three different elongation domains appear: a domain where the elongation is low/medium (A00 5 2.5 i.e. H 5 28 mm) and for which the discrepancy between the two methods is less than looh: the agreement is fair and allows validation of the rheometer measurement: an instability zone where no measurement is possible: this phenomenon is related to the stress-induced crystallization of the elastomer during the biaxial elongation and will be analyzed in Section

5;

in the last domain of large elongation (Aee 2 4.25 i.e. H 2 42 mm), the uncertainty of the magnetic probe measurement increases as a result of local variations of the thickness (the sample slides slightly between the clamps).

6) APPLICATIONS

5.1) Eqdbiaxial Rheological Behavior of an llTRlSBR Elamtomer

We present the equibiaxial stress/strain data for a natural mbber/SBR compound. used by SNECMA-SEP for mechanical applications (flexible pieces) obtained for static inflation at ambient temperature (25°C). This elastomer is a natural rubber (NR) sulfur cured (SEV system) filled with carbon black. For reasons of confidentiality, the exact composition of this material cannot be given. The pressure is imposed at a given value. The different parameters are calculated using the image processing software associated to the bubble inflation rheometer (inflated membrane thickness, curvature radius, hoop stress and stretch). The procedure is repeated until the elastomer limit of elongation is reached (before bursting). We proceed to five measurements for each pressure value. Note that the uncertainty for the measurements is less than 4%, which is quite fair. In Fig. a domain where the elongation is low/medium (A00 5 2.5) and for which the rheological behavior of the elastomer is quasi-linear;

an instability zone where no measurement is possible:

a last domain of large elongation (Aee 2 4.25), where strain hardening occurs.

At small and medium deformations, chain segments of the NR-network are linearly extended. The stressstrain curve is linear. Then, the domain of instability begins and the bubble inflates quickly at constant pressure. Figure 3 illustrates this phenomenon: image 3b shows the bubble just before the instability and image 3c the bubble just after the instability; the elongation goes from 2.5 to 4.25 as the pressure remain constant at 0.6 bar. This phenomenon is well known for inflation of NR membrane and can be attributed to strain-induced crystallization (see Sdon 5.2).

In Fig. 9, the mean curve for equibiaxial stress-elongation of the elastomer studied is plotted and compared with the uniaxial stress-elongation curve. The tensile test was performed with an Instron 1100 equipment using a standard H2 sample 46-002). The elongation speed was fixed to 1 mm/min. In that case, uniaxial and equibiaxial curves are very similar (the maximum discrepancy is 15Oh). This is not the case if we plot the true stress-strain data (Rg. 10). Depending on the type of elongation, the true strain and stress are defined below:

Uniaxial elongation Equibiaxial elongation (9)

We have now for a given value of true stress the equibiaxial true strain that is two times the uniaxial true strain. tests and balloon inflation similar to the phenomenon we have observed. They point out that the instability zone is related to the stress-induced crystallization during elongation at ambient temperature. In addition, they observe that crystallization appears first where stresses and deformations are maximum (the center of the tensile sample and the pole of the bal- loon). Then the crystallization front is progressively propagated in the whole natural rubber. In particular, they point out stress softening for crystallizing elastomers when the Same sample is submitted to successive elongation cycles. This effect is the so-called Mullins effect. We have submitted the natural rubber/ SBR compound elastomer to three successive elongation cycles, and we observe the same behavior as plotted in Rg. I I. This simple experience demonstrates that the instability phenomenon for this elastomer is due to crystallization during elongation. Another effect induced by the crystallization is strain hardening. Because of the stronger alignment of chain segments, crystallization appears preferably at locations of large strain or stress, respectively, and leads to a reinforcement of these weakest domains of the network. The resulting improved tensile strength properties are typical for NR-networks and a consequence of the high cis-portion in this rubber, related to the ability to crystallize, together with the very slow crystal growth velocity (19, 20).

6) CONCLUSION

A bubble-inflation rheometer has been developed and optimized. We have validated the thickness computed using the area change by measuring the thickness using a standard magnetic probe. In addition, we demonstrate that the spherical bubble assumption together with the volume balance assumption will not yield useful values for this application.

Then, biaxial rheological characterization of a natural rubber/SBR compound elastomer was conducted using this rheometer. For this elastomer, the biaxial behavior is very similar to the uniaxial behavior. It should be noted that the comparison depends upon whether we use the true stress and strain or the uniaxial and hoop stress and strain. The stress-induced crystallization phenomenon has been hghlighted.

Future works will allow us to proceed to dynamic inktion by using a high-speed video camera and the relevant image grabber. In order to proceed to polymer inflation, an air convection Oven will be added to the present bubble-inflation rheometer. This will dm us

3 Fig. 2 .

 32 Fig. 2. schematic view of the bubble in_flation rheometer.

524Fig. 3 .

 3 Fig. 3. In-

eFig. 4 .Fig. 6 .

 46 Fig. 4. Video camera images processing.

7 .

 7 , F. M. Schmidt, Y. Le Maoult, M. Rachik, and F. Abbe efeo amparison between magnetic probe and rheometer relative thickness versus bubble height.

Ftg. 8 .Fig. 10 .

 810 Fig. 10. comparison between biavial and uniaxial true stress-strain data for a natuml rubber, SBR compound (T = 25°C).

6

 6 Aee m. I I. Stress-sofining for the elastomer ajter successiue elongation cycles. to heat polymers above the glass transition temperature and to obtain biaxial stress-strain data for polymers. ACKNOWLEDGMENT SEP-Division SNECMA Company supported this research. Special thanks go J. P. kcens, T. Lopez and B. David.

POLYMER ENGINEERING AND SCIENCE, MARCH 2001, Vol. 41, No. 3

531 Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org)


M Rachik

Fabrice Schmidt

N Reuge

Yannick Le Maoult

F Abbe

M Rachik'

F Schmidp

N Reuge3

Y Le Maoull

F Abbe4

Elastomer Biaxial Characterization Using Bubble Inflation Technique. II: Numerical Investigation of Some Constitutive Models

Elastomer biaxial characterization using bubble inflation

technique. II: Numerical investigation of some constitutive models

For a given form of W to be of a practical use, the material parameters have to be determined. In the standard method for material parameter estimations, a specimen is subjected to homogeneous strains (uniaxial, bid extension or pure shear tests) and the stress-strain measurements are used to fit hyper-elastic constants to the measured data. The effectiveness of such a method requires a homogeneous state of deformation during the tests.

In earlier works on non-homogeneous large elastic deformations, Green and Adkins (2) used the Mooney form of W to predict the inflation of a circular plane sheet. The profiles of the inflated membrane were compared with the profiles measured by Treloar (3). This membrane inflation technique was also used by Winman et aL (4) in biomechanics for material identification of soft tissue. The authors proposed an on@nal approach to investigate the domain of definition of W.

During the last decade, many research works were devoted to the material parameter estimation using inverse methods. This approach based on mixed experimental numerical methods does not require the homogeneity of deformation that is needed for the socalled standard method. The inverse identification uses the straight comparison between measured response and the numerical results from the considered constitutive model. This method was successfully used to estimate material parameters of many complex constitutive equations in the inelastic range. In this work, we will investigate some fitting procedures and various strain energy potentials for mechanical behavior characterization of the elastomer studied in part I of this paper. The investigated material is a natural rubber/SBR compound, used by SNECMA-SEP for mechanical applications (flexible parts). The used inverse method is based on the comparison between the measured and the calculated responses for the inflation of a circular membrane. The experimental data are obtained by means of the bubble inflation rheometer previously described.

HYPER-ELASTIC MODEL

In this section, we will briefly recall the fundamental equations of an hyper-elastic model. The mechanical behavior of hyper-elastic materials is completely characterized by means of a strain energy potential that is a function of the deformation tensor invariant: where A,, A, and A, are principal stretches.

Hyper-elastic materials are often nearly or fully incompressible; hence they undergo volume-preserving deformation, which can be expressed as:

z3 = 1 (2)
Among available techniques for taking into account the incompressibility constraint, the Lagrange multiplier is the most convenient one; see Sussman et aL (5) for a complete review on the subject. In its simplest form. the strain energy function of Eq 1 is replaced by: where P is the Lagrange multiplier, equivalent to the hydrostatic pressure.

IDENTIFICATION OF HATERIAL

PARAMETERS

The mechanical characterization of materials is essential in the vast majority of engineering applications. Once the constitutive model is established, it is necessary to determine its material constants. In the so-called direct identification method, the material parameters are determined using measurement on some test specimen with simple boundary condition (uniax- ial extension, torsion etc.). The main drawback of this method is the strain homogeneity requirement, which limits the choice of test specimen and boundary conditions; Michino et al. (6). These difficulties can be overcome using the inverse method for which the strain homogeneity is not required. When used in material parameter identification, the inverse method consists in finding the unknown parameters in such a way that the response calculated with the constitutive model matches the measured response. The method combines optimization techniques with numerical analysis based on finite element or boundary element method to obtain the computed response; Den Camp et aL (7). The inverse technique was successfully used by Aoki et aL (8) to identify constitutive model including damage phenomena, and also by Tillier et aL [START_REF] Tillier | Inverse method for the characterization of mechanical behavior of polymers under biaxial high velocity loading[END_REF] to characterize time dependent constitutive equations. In this section, we compare material constants obtained by direct identification method from uniaxial and equibiaxM extension test data with those obtained by an inverse method from the bubble inflation test data. First, we will briefly recall the standard (direct) identification procedure and then we will describe the inverse method.

1. Direct Identification Method

The mechanical behavior of a rubber-like solid (incompressible) material can be determined by means of experiments involving simple deformation modes like uniaxial extension, equibiaxial extension and pure shear. Once the test data are available, the material constants in the strain energy function can be fitted to these experimental data by means of least square fit procedure. Given 'n" measured stress-&-ain pairs, the best set of constants is the one that miniinizes the following error: Spt: measured nominal stress axisymetric solid elements with hybrid formulation, reduced integration and incompatible modes). Since the obtained results are quite similar, the 2 node axisymetric shell element was adopted for its computa-Sp: theoretical stress derived from a specific strain energy potential

Uniaxial and equibiaxial extension data

For the investigated material, uniaxial and equibiaxial extension test data are shown in Fig. 1.

Depending on the constitutive model, Eq 4 can be non-linear in material parameters. Consequently, a non-linear least squares fit is needed.

Inveme Identification Method

In this work the hyper-elastic constants fitted to uniaxial and biaxial extension test data were used in finite element analysis of the bubble inflation test. We have noticed a substantial difference between the obtained results and the measurements. This discrepancy has motivated the development of an inverse method, which combines optimization technique with finite element analysis of the bubble inflation test described in part I of this paper. tional cost effectiveness. A complete description of the axisymetric finite element model is given in Rg. 2.

The material parameter estimation consists in adjusting the parameters in the finite element model until the calculated displacements u* (bubble height Fig. 3) match the measured displacement il, Gelin et aL (10). It should be pointed out that the comparison between the computed and the measured displacements is performed at the center of the circular membrane. The cost function to minimize with respect to material parameters is given by:

r= rip) = u* - (5)
To reduce the errors on cost function and its gradient, several finite element models were compared ( For physical considerations like Drucker stability, the material parameters must be constrained, Schnur et al. (1 l), which can be written in the form:

C,(p) 50 i= 1, m

(6) These constraints can be taken into account by means of the penalty method, and the modified cost function is then given by:

The cost function @*(p) is minimized by means of the Levenberg-Marquardt algorithm. A detailed description of the used inverse method was given by Pilvin (12).

IDENTIFICATIOlli OF SOME EYPER-ELASTIC MODELS

Since elastomeric materials are increasingly used in industry, great research effort was devoted to the development of constitutive models that describe realis- tically their behavior. Hence several strain energy functions were proposed. In this work, the material parameters of some popular strain energy functions are estimated using the direct identification method and the inverse method. With the direct method, the hyper-elastic constants are fitted to uniaxial and equibiaxial extension test data, whereas with the inverse method, the constants are obtained from the bubble inflation test data. For each constitutive model, the following investigations are performed: Material parameter estimation using direct identification method and uniaxial/equibiaxial extension test data.

Material parameter estimation using inverse method and bubble intlation test data. Check of correlation between the numerical results from the hyper-elastic model and the equibiaxial extension test data.

The finite element calculations are performed using the material parameters obtained by the direct method and those obtained by the inverse method.

w= CI,(Zl -3) + C& -3) (9)
In this work, only the Mooney-Rivlin model (polynomial model with n = 1) is identified. For high order models, oscillations and material instability occur because the available test data are not sufficient.

The material parameter values obtained by direct method, those obtained by the inverse methods, and the associated constraints are given in Table 1.

In Fig. 4, experimental and calculated pressurebubble height curves are represented, along with a comparison between direct identification and inverse identification when the obtained material parameters are used in finite element analysis of the bubble inflation test. It should be noticed that for large strains, the correlation between numerical results from the constitutive model and experimental results is equally poor for the inverse method as for the direct method. But for small strains the inverse method seems to be more effective. The observed discrepancy is due to the low order of the constitutive model (n = 1).

When the material parameters obtained are used in an equibiaxial test simulation, both the inverse method and the direct method give satisfactory correlation between numerical and experimental results (Fis. 3.

Yeoh Model

It has been shown experimentally that the sensitivity of the strain energy potential to changes in the first invariant I, is largely dominant, and the model prediction can be enhanced when neglecting the second invariant contribution that is, in addition, difficult to measure. Following these considerations, Yeoh (13) proposed a reduced form of the polynomial strain energy function where the second invariant contribution is completely omitted:

3 w = 2 c, (I, -3)' (10)
The material parameters obtained by the direct method, those obtained by the inverse methods, and the associated constraints are summarized in Table 1. These material parameters are used in numerical simulation of bubble inflation test. In Flg. 6. experimental and computed pressure-bubble height curves are presented. It is clearly shown that the inverse identification is better than the direct one. We can also notice that the numerical results are in good agreement with experiments for the simulation of the equibiaxial extension test. A comparison between the direct identification and the inverse one is given in Rg. 7. Compared to the polynomial model, this model is more effective because it contains higher powers of the first invariant and the contribution of the second invariant to the strain energy is neglected.

Armda-B~y~e Model

This strain energy function. which depends on the first invariant only, was developed by means of the mechanical statistical theory, Armda et aL (14):

7050'

where p and A, are material parameters to be estimated.

The identified parameter values obtained by the di- rect method, those obtained by the inverse methods, 1 2 3

4 5
The numerical results from this constitutive model are in good agreement with experiments for the simulation of bubble idation test in Fig. 8, as well as for the simulation of equibiardal extension test in Rg. 9.

Once again, the inverse identification method seems to be more effective than the direct one when the obtained material parameters are used in numerical simulation of the bubble inflation test.

Von der Waals Modal

In our work, this model fitted in its general form leads to material instability. Thus a particular form, where the second invariant contribution is omitted, was used:

2 1-33 -(A; -3)(Log( 1 -q) + q) -a( *).> I,-3 f
and the associated constraints are given in Table 1.

where &, p and a are hyper-elastic constants Fig.

6-

5-

4-

3-

Yeoh model

Inverse identlfication

5

The material parameters obtained by the direct method and those obtained by the inverse methods are summarized in Table 1. They are used in numerical simulation of the bubble inflation test in Rg. 10 and the equibiaxial extension test in Rg. 11. For both tests, the numerical results are in good agreement with experiments. With the bubble inflation test, the inverse method is more effective than the direct one, but with the equibiaxial extension test, a discrepancy is observed for large strain. and the associated constraints are summarized in Table 1. When these material parameters are used for numerical simulation of the bubble inflation test in Fig. 12, we can clearly see that the model fails in the membrane behavior prediction when the direct identification is used. The discrepancy becomes significant in the instability zone (27 < bubble height < 43).

For the equibiaxial extension test, the numerical results obtained from the constitutive model are in good agreement with experiments in Rg. 13. 4.6. Ogdem Model 6. CONCLUSION Among the strain energy functions investigated, the Ogden model holds a particular position since the strain energy function is expressed in terms of the principal stretches hl, h2 and A,:

where q and pi are material parameters.

The identified parameter values obtained by the direct method, those obtained by the inverse methods, Experimental data obtained by means of a bubble inflation rheometer are used to characterize the mechanical behavior of a natural rubber. Several hyperelastic constitutive models are considered. Among the strain energy functions investigated in this work, the best fit is obtained when the second invariant contribution is omitted. Except with the Ogden and the Mooney-Rivlin models, the correlation between numerical result from the constitutive model and experiments is satisfactory for the bubble inflation test as well as for the equibiaxial extension test. In addition, it can be noticed that the model prediction is enhanced when the second invariant contribution is neglected. An inverse identification method, which combines optimization technique and finite element analysis, is used to estimate material parameters. The parameter values obtained are compared with those obtained by the standard direct method. With regard to the error analysis, we observed the same tendency for the various constitutive models. To illustrate this, we present the error in numerical results from the simulation of the bubble inflation test and those from the equibiaxial extension test simulation with the Yeoh model. In particular concerning the bubble inflation test, the evolution of the error in the pressure with respect to the bubble height, for both the direct identification and the inverse one, is given in F@. 14. The material parameters obtained by the direct identification seem to fail in the description of non-homogeneous state of deformation.

The same investigations were carried out for the equibiaxial extension test. As shown in Fig. 15, presenting the error in nominal stress with respect to the nominal strain, the comparison between the two methods shows that the error remains limited in the case of the direct identification, whereas it tends to grow in the case of the inverse identification.

In light of the previous observations, the inverse method seems to be a good tool for the material parameter estimation. Although such a method does not need the restrictive requirement of the strain state homogeneity, it can be further improved or calibrated when associated with simple homogeneous test data, like extension or equibhxial extension.

( 4 ) 3 Fig. 1 .

 431 Fig. 1. Nominal strain-nominal stress mmes obtainedfrom uniaxial and equibiuxial extension tests.

  f[, , , , , height(mm) , Check of correlation between the numerical results from the hyper-elastic model and bubble inflation test data.

4. 1

 1 Polynomial Strain Enew Function and 1-Variant This model is based on phenomenological theory first initiated by Rivlin, who proposed a polynomial function as the general form of the strain energy potentiak n where Cij are hyper-elastic constants. This general form includes some popular energy hctions like the Mooney-Rivlin model and the Neo-Hookean one. The Mooney-Rivlin function is obtained by setting n = 1 and C, = 0 except C,, and C,,:

Rg. 4 .

 4 Comparison of numerical results from Mooney-Rslin constitutive model (direct and inverse identijications) and bubble in@tion test data.

3 Rg. 5 .

 35 Comparison of numerical results from Mooney-Rivlin constitutive model (direct and inverseidentifications) and equibiaxial atension test data.

5 flg. 10 .

 510 Comparison of numerical resulkfrom Van der Waals constitutive model [direct and inverse identijications) and bubble injlation test data.

Fig. 11 .

 11 Fig. 11. Comparisonof numerical resultsfrom Van der Waals constitutive model (direct and inverse identiiaths) and equibiaxid extension test data

3 Q. 15 .

 315 Error on numerical results of equibiaxiat test simulation with Yeoh model (comparison between direct and inverse identiiations).

3 Fig. 3. Pressure-bubble height twn rheometer and used for

2 <p> = <pl p2....pn> material parameters and 3 node axisymetric shell elements, 4 and 8 node N, is the number of measured bubble heights. FE model 534 FYg. 2. Awisymetricjinite element model of bubble i n ~n test used with inverse identiiation method POLYMER ENGINEERING AND SCIENCE, MARCH 2001, Vol. 41, No. inverse identiiation. c~uve Obtained by the bubble hjkt 0 9 is a weighting matrix used to improve the result accuracy when more than one test data are available (bubble height, thickness.. .)

Table . 1. Material Parameter Values for Different Models.

 .

	Model	Parameter	Direct identification	Inverse identification	Constraints
	polynomial (n = 1)	CIO (MP4 COl (MP4	0.101732 1.243427 X lo-*	0.1 571 831 1.010127 X 10"	CIO + COl > 0
	Yeoh	ClO (MP4 cm (Mpa) c , (MP4	0.1 371 83 1.415 X lo9 1.5568 X	0.19405 0.93719 X lW5 0.33713 X lo4	ClO > 0 c20 > 0 c,>o
	Arruda-Boyce		0.281 91 4.695	0.327993 4.4	P>O A, > 0
	Van der waals		0.244314	0.384766
			11.2545 -7.9445 x 10-02	10.25569 0.189963
	Ogden (n = 2)	P1 (MP4 P2 (MP4	0.219472 1.79949 x 104	0.270533 0.131425 X	P1 + P2 > 0
		a1	-0.7487	-0.7266249
		a2	-2.82844	-3.588489

6. Comparison of numerical results from Yeoh constitutive model (direct and inverse identii- cations) and bubble inflation test data

	8-
	7-
	Fig. 7. Comparison of numerical
	results from Yeoh constitutive
	model (direct and inverse identiji-
	cations) and equibiawiaI extension
	test data.
	Fig. 8. Comparison of numerical
	results from Anuda-Boyce consti-
	tutive model (direct and inverse
	identifications) and bubble inta-
	tion test data
								(12f"l
	Fig. 9. Comparison of numerical results from Arruda-Boyce consti-			Yeoh model
	tutive model (direct and inverse
	identiiationss) and equibiawial ex-
	tension test data
				-------
	Direct identification 0,m y , I ,o ExQeriment;	height(mm)
	0
	0	10	20	30	40	50	60	70	80

POLYMER ENGINEERING AND SCIENCE, MARCH 2001, Vol. 41, No. 3

	0,l I
	FQ. 12. Comparison of numerical
	resultsfrorn @den In = 2) consti-
	tutive model (direct and inverse identilications) and bubble inm-tion test data.
	FYg. 13. Comparison of Rumerical
	resultsfrorn Ogden (n = 2) consti-
	tutive model (direct and inverse
	identtfrcations) and equibiaxial ex-
	tension test data
	FQ. 14. Error in numerical results
	of bubble infrcltion test simulation
	with Yeoh model (comparison be-
	tween direct and inverse identit-
	mtionsl.

541 Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org)


Ali Feizmohammadi

Joonas Ilmavirta
email: joonas.ilmavirta@jyu.fi

AND Yavar Kian
email: yavar.kian@univ-amu.fr

Lauri Oksanen
email: l.oksanen@ucl.ac.uk

RECOVERY OF TIME DEPENDENT COEFFICIENTS FROM BOUNDARY DATA FOR HYPERBOLIC EQUATIONS

Keywords: Dirichlet to Neumann map, Gaussian beam, inverse problems, lens data, light ray transform, magnetic potential

We study uniqueness of the recovery of a time-dependent magnetic vector-valued potential and an electric scalar-valued potential on a Riemannian manifold from the knowledge of the Dirichlet to Neumann map of a hyperbolic equation. The Cauchy data is observed on time-like parts of the space-time boundary and uniqueness is proved up to the natural gauge for the problem. The proof is based on Gaussian beams and inversion of the light ray transform on Lorentzian manifolds under the assumptions that the Lorentzian manifold is a product of a Riemannian manifold with a time interval and that the geodesic ray transform is invertible on the Riemannian manifold.

1. Introduction 1.1. Statement of the problem. Let (M, ḡ) be a 1 + n dimensional Lorentzian manifold with boundary. Throughout this paper, we will assume that (M, ḡ) has a global product structure, that is to say M = [0, T ] × M , ḡ = -dt 2 + g, where T > 0 and (M, g) denotes a smooth compact connected Riemannian manifold of dimension n 2, with smooth boundary ∂M . We assume that g ∈ C 6 (M ; Sym 2 M ), where Sym 2 M denotes the bundle of symmetric two-tensors over M . We denote by ∆ ḡ the Laplace-Beltrami operator given by ∆ ḡ := div ḡ ∇ ḡ, where div ḡ (resp., ∇ ḡ) denotes the divergence (resp., gradient) operator on (M, ḡ). In local coordinates (x 0 := t, x 1 , . . . , x n ) and for each u ∈ C 2 (M) we have:

∆ ḡu = n i,j=0 |ḡ| -1/2 ∂ x i (|ḡ| 1/2 ḡij ∂ x j u),
where ḡ-1 := (ḡ ij ) 0 i,j n and |ḡ| := |det ḡ|. As ḡ = -dt 2 +g, we have ∆ ḡ = -∂ 2 t +∆ g , where ∆ g is defined analogously. Consider a complex-valued scalar function q(t, x) (electric potential) and a complex-valued one-form A (magnetic potential), such that the following regularity assumptions are satisfied q ∈ C(M) and A ∈ C 1 (M; T * M).

(1.1)

In local coordinates, the one-form A can be expressed as

A(t, x) = b(t, x) dt + n i=1 ω j (t, x) dx j = b(t, x) dt + ω(t, x), (1.2)
where ω is a time-dependent one-form on (M, g). Given A and q as above, We consider the initial boundary value problem (IBVP)

   -∆ ḡu + A∇ ḡu + qu = 0, on M, u = h, on (0, T ) × ∂M, u(0, •) = 0, ∂ t u(0, •) = 0 on M, (1.3) with non-homogeneous Dirichlet data h ∈ H 1 0 ((0, T ) × ∂M ). Note that in local coordinates, A∇ ḡu = -b ∂ t u + n i,j=1 g ij ω i ∂ j u. We introduce ν, the outward unit normal vector to (0, T ) × ∂M , and define the hyperbolic Dirichlet-to-Neumann (DN in short) map, given by

Λ A,q : h → (∂ ν u - (Aν)u 2 ) |(0,T )×∂M
where u solves problem (1.3). This is well-defined as equation (1.3) admits a unique solution u ∈ C(0, T ; H 1 (M )) ∩ C 1 (0, T ; L 2 (M )), with ∂ ν u| (0,T )×∂M ∈ L 2 ((0, T ) × ∂M ). This follows from [29, Theorem 2.1] as explained in Section 2.1. The goal of this paper is to study the unique recovery of the complex valued coefficients A and q given Λ A,q , up to the natural obstructions discussed in the next section.

1.2. Natural obstructions. The first obstruction concerns the recovery of the magnetic potential A. Indeed, for j = 1, 2, fix (A j , q j ) defined as above and assume that there exists ψ(t, x) ∈ C 2 (M) with ψ| (0,T )×∂M = 0, such that

A 1 = A 2 + 2 dψ, q 1 = q 2 + ∆ ḡψ -A 2 ∇ ḡψ -∇ ḡψ, ∇ ḡψ ḡ , (1.4)
where d denotes the exterior derivative on M, that acts on f ∈ C ∞ (M) through df = ∂ t f dt + df with d denoting the exterior derivative on M , and •, • ḡ denotes the inner product on (M, ḡ). Then, for u j the solution of (1.3) with A = A j , q = q j , j = 1, 2, it holds that u 1 = e ψ u 2 and, using the fact that ψ |(0,T )×∂M = 0, we obtain

∂ ν u 1 - (A 1 ν)u 1 2 = ∂ ν u 2 + (∂ ν ψ)u 2 - (A 1 ν)u 2 2 = ∂ ν u 2 - (A 2 ν)u 2 2 .
This proves that Λ A 1 ,q 1 = Λ A 2 ,q 2 , but A 1 = A 2 as soon as ψ does not vanish identically. In other words, the DN map Λ A,q is invariant with respect to the gauge transformation given by (1.4) and the best we can expect is the recovery of the coefficients A and q from Λ A,q modulo the gauge invariance (1.4).

22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE

The second obstruction to our problem is due to finite speed of propagation for the wave equation. Let us define the set

D = {(t, x) ∈ M | dist(x, ∂M ) < t < T -dist(x, ∂M )}.
Due to domain of dependence arguments (see [START_REF] Kian | Unique determination of a time-dependent potential for wave equations from partial data[END_REF]Section 1.1] and [26, Section 1.1]), it is not possible to recover the restriction of any of the coefficients A and q on the set M \ D from Λ A,q . Therefore, for our problem, the best we can expect, is to recover the coefficients modulo the gauge invariance above, on the set D.

1.3. Main result. Before stating the main result, we need to recall the definition of the geodesic ray transform on the (spatial) Riemannian manifold M . For each (y, v) ∈ SM , with SM denoting the unit tangent bundle of M , let γ(•; y, v) denote the unit speed geodesic starting at point y, in the direction v, that is: ∇ g γ γ(•; y, v) = 0, γ(0; y, v) = y, γ(0; y, v) = v. For any (y, v) ∈ SM , we define τ exit (y, v) through: τ exit (y, v) := inf {t > 0 | γ(t; y, v) ∈ ∂M, γ(t; y, v) / ∈ T γ(t;y,v) ∂M }.

We now define

∂ -SM := {(y, v) ∈ SM | y ∈ ∂M, v, ν(y) g < 0, τ exit (y, v) < ∞}, (1.5)
where ν denotes the outward normal unit vector on ∂M . Henceforth, for the sake of brevity, we use the term maximal geodesic to refer to the geodesics γ(•; y, v) (or γ(•) in short) with (y, v) ∈ ∂ -SM , over their maximal interval of definition in M int , that is the interval I := (0, τ exit ).

Definition 1.1. Let (y, v) ∈ ∂ -SM and let γ(•; y, v) : I → M . We define the geodesic ray transform of (f, α) ∈ C(M ) × C(M ; T * M ), as follows:

I γ (f, α) := I [f (γ(t)) + α(γ(t)) γ(t)] dt.
We also need to recall the definition of the solenoidal component, α s , of a oneform α with local representation α

= n k=1 α k dx k . Let δα := n i,j=1 1
√ g ∂ i ( √ gg ij α j )
denote the divergence operator on M sending one-forms to functions. Given any α ∈ L 2 (M ; T * M ), it can be uniquely decomposed as

α = α s + dψ, (1.6)
where δα s = 0 and ψ ∈ H 1 0 (M ) solves ∆ g ψ = δα in the weak sense on M . This is called the Helmholtz decomposition (see for example [START_REF] Taylor | Partial Differential Equations[END_REF]). We will be working with C 1 (M ) one-forms. In this case, one can immediately see that since δα ∈ C(M ) ⊂ L n (M ), elliptic regularity implies that ψ ∈ W 2,n (M ) ⊂ C 1 (M ). We will use this 22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE observation later in the paper to derive the smoothness properties for the gauge. With these notations, we can now state the main geometric assumption on (M, g). Hypothesis 1.2. Let (M, g) be a compact connected Riemannian manifold with smooth boundary. We say that the geodesic ray transform is injective on M with respect to functions f ∈ C(M ) and one-forms α ∈ C(M ; T * M ), if the following holds:

I γ(•;y,v) (f, α) = 0, ∀(y, v) ∈ ∂ -SM implies that f ≡ 0 and α s ≡ 0.
According to Theorems 3 and 4 in [START_REF] Stefanov | Stability estimates for the X-ray transform of tensor fields and boundary rigidity[END_REF], Hypothesis 1.2 will be fulfilled if M is simple. This condition can also be fulfilled by a non-simple manifold. We refer to Section 2.2 for a more detailed discussion about this aspect.

Finally, let us introduce the set E ⊂ D where we recover the coefficients. For any x ∈ M , we define

D g (x) := sup {τ exit (y, v) | (y, v) ∈ ∂ -SM, x is in γ(•; y, v)} and let D g (M ) := sup {D g (x) | x ∈ M }. For T > 2D g (M ), we define E := {(t, x) ∈ M | D g (x) < t < T -D g (x)}. Theorem 1.3. Let g ∈ C 6 (M ; Sym 2 M ), A 1 , A 2 ∈ C 1 (M; T * M) and q 1 , q 2 ∈ C(M). Assume that supp (A 1 -A 2 ) ⊂ E, supp (q 1 -q 2 ) ⊂ E, and A 1 (t, x) = A 2 (t, x), ∀(t, x) ∈ (0, T ) × ∂M.
If Hypothesis 1.2 holds, then Λ A 1 ,q 1 = Λ A 2 ,q 2 implies that there exists ψ ∈ C 2 (M) with ψ| (0,T )×∂M = 0, such that (1.4) holds.

The proof of this theorem relies in part on injectivity of the so-called light ray transforms of vector valued functions. Recall that a curve β in M is a null geodesic (or a light ray), if ∇ ḡβ β = 0 and β, β ḡ = 0. Given the product structure of M, we can parametrize maximal null geodesics through β(t) = ( t + t, γ(t)), t ∈ I where γ is a unit speed maximal geodesic in M and t ∈ R. We define the light ray transforms L q and L A of q ∈ C(M) and A ∈ C(M; T * M) as follows:

L q (β) := I q(β(t)) dt and L A (β) := β A.
We have the following proposition, that will be proved in Section 5.

Proposition 1.4. Let (q, A) ∈ C(M) × C 1 (M; T * M) be such that supp q, supp A ⊂ E. If Hypothesis 1.2 holds, then we have the following statements:

(i) If L q (β) = 0 for all maximal null geodesics β in M, then q ≡ 0 on M.

(ii) If L A (β) = 0 for all maximal null geodesics β in M, then there exists ψ ∈ C 2 (M) vanishing on ∂M such that A ≡ dψ on M.

1.4. Previous literature. Results related to the recovery of coefficients for hyperbolic equations can in general be divided into two categories of time-independent and time-dependent coefficients. Starting with the seminal works [START_REF] Belishev | An approach to multidimensional inverse problems for the wave equation[END_REF][START_REF] Belishev | To the reconstruction of a Riemannian manifold via its spectral data (BC-method)[END_REF], there is an extensive literature related to the recovery of time-independent coefficients for hyperbolic equations. These results usually rely on the Boundary Control method, developed in [START_REF] Belishev | An approach to multidimensional inverse problems for the wave equation[END_REF][START_REF] Belishev | To the reconstruction of a Riemannian manifold via its spectral data (BC-method)[END_REF] and a time sharp unique continuation theorem [START_REF] Tataru | Unique continuation for solutions to PDE; between Hörmanders theorem and Holmgrens theorem[END_REF], which provide the building blocks of very general results. We refer the reader to [START_REF] Kian | Application of the boundary control method to partial data Borg-Levinson inverse spectral problem[END_REF] for an introduction to the method and to [START_REF] Kurylev | Inverse problems for the connection Laplacian[END_REF] for an example of a state of the art result in this direction. We also refer to [START_REF] Belishev | Recent progress in the boundary control method[END_REF][START_REF] Katchalov | Inverse boundary spectral problems[END_REF] for review. The unique continuation theorem in [START_REF] Tataru | Unique continuation for solutions to PDE; between Hörmanders theorem and Holmgrens theorem[END_REF] fails if the dependence of the coefficients on time is non-analytic and therefore extension of these results for general time dependent coefficients is not possible (see e.g. [START_REF] Alinhac | Non-unicité du probléme de Cauchy[END_REF][START_REF] Alinhac | A non uniqueness result for operators of principal type[END_REF]). We refer the reader to [START_REF] Eskin | Inverse hyperbolic problems with time-dependent coefficients[END_REF] for a uniqueness result, when the dependence of the coefficients on time is real analytic. Starting with [START_REF] Bukhgeim | Global uniqueness of a class of multidimensional inverse problem[END_REF], methods based on Carleman estimates have also been quite fruitful in deriving uniqueness results for time-independent coefficients of hyperbolic equations. Contrary to the Boundary Control method, where the best known stability estimates are double logarithmic [START_REF] Bosi | Stability of the unique continuation for the wave operator via Tataru inequality and applications[END_REF], these methods tend to give strong stability estimates. We also mention [START_REF] Liu | A Lipschitz stable reconstruction formula for the inverse problem for the wave equation[END_REF] where Boundary Control method is combined with complex geometric optics and stronger estimates are obtained for low frequencies.

In the time-dependent category, most of the results are obtained for the reconstruction of the zeroth order term, q, and are based on a use of geometric optic solutions for the wave equation. Let us mention that this approach has also been used in the time-independent category [START_REF] Bellassoued | Stability estimates for the anisotropic wave equation from the Dirichlet-to-Neumann map[END_REF][START_REF] Bellassoued | Lipschitz stability for a hyperbolic inverse problem by finite local boundary data[END_REF][START_REF] Kian | Stability of the determination of a coefficient for wave equations in an infinite waveguide[END_REF][START_REF] Stefanov | Stable determination of the hyperbolic Dirichlet-to-Neumann map for generic simple metrics[END_REF] to obtain strong stability estimates although they suffer from considerably stronger geometric assumptions compared to the Boundary Control method. Methods based on geometric optics were used in the context of recovery of time-dependent coefficients starting with [START_REF] Stefanov | Uniqueness of the multi-dimensional inverse scattering problem for time dependent potentials[END_REF]. Among the literature of results in this direction, we refer to [START_REF] Aicha | Stability estimate for hyperbolic inverse problem with time-dependent coefficient[END_REF][START_REF] Isakov | Completness of products of solutions and some inverse problems for PDE[END_REF][START_REF] Kian | Stability in the determination of a time-dependent coefficient for wave equations from partial data[END_REF][START_REF] Ramm | An inverse problem of the wave equation[END_REF][START_REF] Salazar | Determination of time-dependent coefficients for a hyperbolic inverse problem[END_REF]. The leading coefficients for the wave equation in all these results are constant. Uniqueness of zeroth order coefficient q for a variable coefficient wave equation was considered in [START_REF] Kian | Recovery of time-dependent coefficient on Riemanian manifold for hyperbolic equations[END_REF] where the recovery of the potential was based on inversion of geodesic ray transform for scalar functions. It should be noted that even in the case where A = 0, the result in this paper is a significant improvement of [START_REF] Kian | Recovery of time-dependent coefficient on Riemanian manifold for hyperbolic equations[END_REF], since there (M, g) was assumed to be simple.

Approaches based on global geometric optic solutions fail if the Riemannian manifold (M, g) is not simple (see Section 3.5). This motivates the use of Gaussian beams in the current paper. Gaussian beams were introduced in [START_REF] Babich | The complex space-time ray method and quasi-photons[END_REF][START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF] and they were first used in the context of inverse problems in [START_REF] Belishev | Boundary control and quasi-photons in the problem of a Riemannian manifold reconstruction via its dynamical data[END_REF][START_REF] Katchalov | Multidimensional inverse problem with incomplete boundary spectral data[END_REF]. We refer the reader to [START_REF] Katchalov | Inverse boundary spectral problems[END_REF] for a thorough presentation in the case of a wave equation with a smooth metric, a smooth electric potential and no magnetic potential. This paper is concerned with the reconstruction of time-dependent vector valued coefficients for the 22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE wave equation under weak geometrical assumptions on the spatial manifold (M, g) and weaker regularity assumptions on the coefficients. We use Gaussian beams to reduce the inverse problem to the inversion of the light ray transform of the unknown coefficients. The closest previous work to this reduction is [START_REF] Stefanov | The inverse problem for the Dirichlet-to-Neumann map on Lorentzian manifolds[END_REF], where the authors study the problem of recovery of the geometry along with a time-dependent magnetic potential A and an electric potential q in a Lorentzian manifold from a micro-local formulation of a Cauchy data set on the boundary. It is shown that if ḡ, A, q belong to some C k , with k sufficiently large, then this Cauchy data set uniquely determines the scattering relation of ḡ along with light ray transforms of A, q. Their approach is based on the study of Fourier Integral Operators and propagation of singularities. The inversion of the light ray transform on a general Lorentzian manifold is left as an open problem. Our Gaussian beam construction makes the reduction to the light ray transform more explicit in terms of the smoothness required. We also succeed in the inversion of the light lay transform, in the sense of Proposition 1.4. Our inversion method for the light ray transform was inspired in part by techniques developed in the context of the Calderón problem [START_REF] Cekić | The Calderón problem for connections[END_REF].

1.5. Outline of the paper. This paper is organized as follows. In Section 2, we discuss the forward problem (1.3) and also discuss the Hypothesis 1.2. In Section 3, we present the Gaussian beam construction near a null geodesic in M. In Section 4, we show the reduction step from the knowledge of the Dirichlet to Neumann map Λ A,q , to the knowledge of the light ray transforms of A, q and conclude that Theorem 1.3 follows from Proposition 1.4. Finally, Section 5 is concerned with the proof of Proposition 1.4.

Preliminaries

2.1. Direct problem. Let Xu := A∇ ḡu + qu, where A and q satisfy (1.1). We consider the wave equation

     -∆ ḡu + Xu = F, in M, u| x∈∂M = h, on (0, T ) × ∂M , u| t=0 = u 0 , ∂ t u| t=0 = u 1 , on M (2.1) It is classical that u is in the energy space C([0, T ]; H 1 (M )) ∩ C 1 ([0, T ]; L 2 (M )) (2.2) when h = 0, F ∈ L 2 (M), u 0 ∈ H 1 0 (M ) and u 1 ∈ L 2 (M ). The wave equation      -∆ ḡv = F, in M, v| x∈∂M = h, on (0, T ) × ∂M , v| t=0 = u 0 , ∂ t v| t=0 = u 1 , on M , (2.3)
22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE was considered in [START_REF] Lasiecka | Non homogeneous boundary value problems for second order hyperbolic operators[END_REF]. It was shown there that if F and u 1 are as above, and u 0 ∈ H 1 (M ) and h ∈ H 1 ((0, T ) × ∂M ) satisfy the compatibility condition

h| t=0 = u 0 | x∈∂M , (2.4)
then the solution v is the energy space (2.2), and

∂ ν v| x∈∂M ∈ L 2 ((0, T ) × ∂M ).
Let us now set u = v -w where v is the solution of (2.3) with F , u 0 , u 1 and h as above, and w is the solution of (2.1) with F = Xv ∈ L 2 (M), u 0 = 0, u 1 = 0 and h = 0. Then u satisfies (2.1) with the same F , u 0 , u 1 and h as in (2.3) for v. As both v and w are in (2.2), so is u. But then

-∆ ḡu = F -Xu ∈ L 2 (M) and ∂ ν u| x∈∂M ∈ L 2 ((0, T ) × ∂M ).
It is straightforward to turn this regularity result to the corresponding estimate

u C([0,T ];H 1 0 (M ))∩C 1 ([0,T ];L 2 (M )) + ∂ ν u L 2 ((0,T )×∂M ) (2.5) ≤ C( F L 2 (M) + h H 1 ((0,T )×∂M ) + u 0 H 1 (M ) + u 1 L 2 (M ) ),
for solutions u of (2.1) under the compatibility condition (2.4).

We write If (t) = t 0 f (s)ds, and show now that u satisfies the estimate

u L 2 (M) ≤ C IF L 2 (M) , (2.6)
when u 0 , u 1 and h vanish identically. The map taking F to u is continuous from L 2 (M) to H 1 (M), and by considering its adjoint, we obtain also continuity from

H -1 (M) to L 2 (M). Let F ∈ L 2 (M)
and define v and w as above, but with u 0 , u 1 and h vanishing also for v.

Then -∆ ḡI v = IF and v L 2 (M) = ∂ t Iv L 2 (M) ≤ C IF L 2 (M) .
Moreover,

w L 2 (M) ≤ C Xv H -1 (M) ≤ C v L 2 (M) ,
and the above two estimates imply (2.6).

2.2.

On inversion of the geodesic ray transform. We are not aware of results on inversion of the geodesic ray transform I γ (f, α) when the metric tensor g is assumed to be only C 6 -smooth. Therefore, for the purposes of the present section, we assume g to be C ∞ -smooth.

It is well-known that simple manifolds satisfy Hypothesis 1.2, see Theorems 3 and 4 in [START_REF] Stefanov | Stability estimates for the X-ray transform of tensor fields and boundary rigidity[END_REF]. It is also likely that the method to invert the geodesic ray transform using convex foliations, originating from [START_REF] Uhlmann | The inverse problem for the local geodesic ray transform[END_REF], can be used to show that Hypothesis 1.2 holds under the assumptions that the boundary of M is strictly convex in the sense of the second fundamental form, and that there is a strictly convex function on M . In [START_REF] Paternain | The geodesic X-ray transform with matrix weights[END_REF] this is shown under the further assumption that f and α are smooth. Let us point out that, even when A = 0, combining Theorem 1.3 and [START_REF] Paternain | The geodesic X-ray transform with matrix weights[END_REF] gives a result on 22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE unique determination of smooth q that falls outside the scope of the closest previous results [START_REF] Kian | Recovery of time-dependent coefficient on Riemanian manifold for hyperbolic equations[END_REF].

Let us now describe a non-simple case satisfying Hypothesis 1.2 as studied in [START_REF] Stefanov | Integral geometry on tensor fields on a class of non-simple Riemannian manifolds[END_REF]. There it is assumed that (M, g) satisfies the following:

(i) M and ∂M have real analytic atlases.

(ii) There is an open set of simple geodesics Γ on a slightly larger manifold ( M , g)

such that T * M ⊂ {N * γ | γ ∈ Γ}. (iii) Any path in M connecting two boundary points is homotopic to a path of the form c 1 ∪ γ 1 ∪ c 2 ∪ γ 2 ∪ • • • ∪ γ k ∪ c k+1
, where c j are paths on ∂M and γ j = γj | M for some γj ∈ Γ. Moreover, γ j intersects ∂M only at its endpoints and is transversal to ∂M . (iv) g is generic in the sense of [START_REF] Stefanov | Integral geometry on tensor fields on a class of non-simple Riemannian manifolds[END_REF]Corollary 1].

Here γ ∈ Γ being simple means that the endpoints of γ are in M \ M and there are no conjugate points on γ, simple geodesics are given topology in the sense of the parametrization (2) in [START_REF] Stefanov | Integral geometry on tensor fields on a class of non-simple Riemannian manifolds[END_REF], and N * γ denotes the conormal bundle of γ, viewed as a 1-dimensional submanifold of M .

Gaussian Beam Solutions

The goal of this section is to construct the so called Gaussian beam solutions

u j ∈ C([0, T ]; H 1 (M )) ∩ C 1 ([0, T ]; L 2 (M )) for j = 1, 2 of the problems -∆ ḡu 1 + A 1 ∇ ḡu 1 + q 1 u 1 = 0, (t, x) ∈ M, u 1 (0, x) = ∂ t u 1 (0, x) = 0, x ∈ M, -∆ ḡu 2 -A 2 ∇ ḡu 2 + (-δA 2 + q 2 )u 2 = 0, (t, x) ∈ M, u 2 (T, x) = ∂ t u 2 (T, x) = 0, x ∈ M, (3.1) taking the form u 1 (t, x) = e iρϕ(t,x) v 1 (t, x) + R 1,ρ (t, x), (t, x) ∈ M, (3.2)
u 2 (t, x) = e -iρ φ(t,x) v2 (t, x) + R 2,ρ (t, x), (t, x) ∈ M, (3.3) with ρ > 1.
Here, δ denotes the divergence operator on (M, ḡ) sending one-forms to functions. The two equations in (3.1) are in essence formal adjoints of each other, with respect to the real L 2 (M) inner product. The phase function ϕ will be chosen so that both the oscillatory parts e iρϕ and e -iρ φ remain bounded in L 2 (M) as ρ → ∞ and such that the principal terms |e iρϕ(t,x) v k (t, x)| are concentrated near a fixed maximal null geodesic in D. The remainder terms R 1,ρ , R 2,ρ will vanish in the limit ρ → ∞.

22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE 3.1. Fermi Coordinates. We will start by reviewing Fermi coordinates near a fixed maximal null geodesic β : [τ -, τ + ] → M, where we are using the time coordinate t as the parametrization for the null geodesic. Here, β(τ -), β(τ + ) denote the start and end points of the maximal null geodesic on the boundary (0, T ) × ∂M . Note, in particular, that β(t) ∈ D for all t ∈ [τ -, τ + ]. Fermi coordinates were first introduced by E. Fermi [18]. In this paper, the geometry has a product structure which makes the construction of Fermi coordinates slightly easier. This is to some extent similar to [START_REF] Santos Ferreira | The Calderón problem in transversally anisotropic geometries[END_REF][START_REF] Kenig | The Calderón problem with partial data on manifolds and applications[END_REF], where a coordinate construction was carried out in the context of an elliptic partial differential equation on a Riemannian manifold with a product structure. We will therefore follow [START_REF] Santos Ferreira | The Calderón problem in transversally anisotropic geometries[END_REF] with some modifications.

Let us introduce notation that will be fixed throughout the remainder of this paper. We begin by embedding (M, g) into a closed manifold ( M , g) and extend the null geodesic β to M := (0, T ) × M such that β(t) is well-defined on the interval [τ --ε, τ + +ε] with ε > 0 a small constant. We also consider extensions of the metric g and the coefficients A 1 , A 2 , q 1 , q 2 to the bigger set M such that the extended metric g is C 6 smooth and the extensions of the coefficients satisfy the regularity assumptions in (1.1) with M replaced by M. Finally, for the sake of convenience, we define the constants a, b, a 0 , b 0 , s -and s + as follows

a = √ 2(τ --ε), b = √ 2(τ + + ε), a 0 = √ 2(τ -- ε 2 ), b 0 = √ 2(τ + + ε 2 ), (3.4)
and

s -= √ 2τ -, s + = √ 2τ + .
We will now present Fermi coordinates near the null geodesic β in M. In all the following arguments, β(t) denotes the parametrization of the null geodesic with respect to the time coordinate in M.

Lemma 3.1. (Fermi coordinates) Let β : (τ --ε, τ + + ε) → M be a null geodesic as above. There exists a coordinate neighborhood

(U, Φ) of β([τ --ε 2 , τ + + ε 2 ]) denoted by (z 0 := s, z 1 := r, z 2 , . . . , z n ) such that: • Φ(U ) = (a, b) × B(0, δ )
where a, b are given by (3.4) and B(0, δ ) denotes a ball in R n with a sufficiently small radius δ only depending on the geometry (M, ḡ) and ε.

• Φ(β(t)) = ( √ 2t, 0, . . . , 0 n times ) for all t ∈ (τ --ε, τ + + ε).
Furthermore the metric ḡ in this coordinate system is C 4 smooth and satisfies

ḡ| β = 2dsdr + n j=2
(dz j ) 2 and ∂ḡ jk ∂z i | β = 0 for all 0 i, j, k n.

22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE

Proof. Let us begin by defining Π : M → M through Π(t, x) = x. Note that Πβ = γ where γ is a unit speed geodesic passing through the point x 0 = Π(β(a)). We choose {α 2 , ..., α n } such that the set { γ(x 0 ), α 2 , ..., α n } forms an orthonormal basis for T x 0 M . Let y 1 denote the arc length parameter along the geodesic γ from the point x 0 . For each 2 k n, let e k (y 1 ) ∈ T γ(y 1 ) M denote the parallel transport of α k along γ to the point γ(y 1 ). Since γ is also parallel along γ, the set { γ(y 1 ), e 2 (y 1 ), ..., e n (y 1 )} forms an orthonormal basis for T γ(y 1 ) M . We now define the coordinate system (y 0 , . . . , y n ) through F 1 : R n+1 → M:

F 1 (y 0 := t, y 1 , ..., y n ) = (t, exp γ(y 1 ) ( n α=2
y α e α (y 1 ))),

where exp p (•) denotes the exponential map on M at the point p. Let us remark that since g ∈ C 6 (M ; Sym2 M ), the map F 1 is locally in C 5 (see for example [START_REF] Deturck | Some regularity theorems in Riemannian geometry[END_REF]). We now define (s := z 0 , r

:= z 1 , . . . , z n ) = F 2 (y 0 , . . . , y n ) through    s := z 0 := 1 √ √ 2 (-t + y 1 ) + a 2 , z j := y j ∀j 2. (3.5)
For the sake of brevity, we will also use the notations z = (s, z ) = (s, r, z ) for this coordinate system. Let us consider the composition map F : R 1+n → M given by

F = F 1 • F -1 2 . It is clear that for all t ∈ [τ --ε, τ + + ε]: F( √ 2t, 0) = F 1 (t, t -a √ 2 2 , 0) = F 1 (t, t -(τ --ε), 0) = β(t), implying that F(s, 0) is injective for all s ∈ (a, b) as β is not self-intersecting on the time interval [τ --ε, τ + + ε]. Furthermore, for all s ∈ (a, b) it holds that ∂ ∂s F(s, 0) = 1 √ 2 (∂ t + γ( √ 2 2 (s -a))), ∂ ∂r F(s, 0) = 1 √ 2 (-∂ t + γ( √ 2 2 (s -a))), ∂ ∂τ F(s, τ v α )| τ =0 = e α ( √ 2 2 (s -a)),
where v α denotes the αth coordinate vector in R n-2 . Thus F(s, z ) is a locally C 5 map in a neighborhood of the null geodesic β such that F(s, 0) is injective and DF(s, 0) is invertible. The inverse mapping theorem applies to deduce that F(s, z ) is a diffeomorphism on a neighborhood of [a 0 , b 0 ] × B(0, δ ) with δ sufficiently small. We then choose Φ = F -1 . Note that Φ ∈ C 5 near the null geodesic, and since Let us now study the form of the metric in this single coordinate chart (U, Φ) given by the z-coordinates. We will first derive the form of the metric in y-coordinates which is just an affine transformation of the z coordinates (the linear part of this affine transformation is unitary). To find the form of the metric in y-coordninates, we note that F 1 preseves the product structure on M and therefore it suffices to check the form of the Riemannian metric g near the geodesic γ = Π(β) in M . Let the indices i, j, k run between 1 and n and the indices α, β between 2 and n. Since the set { γ(y 1 ), e 2 (y 1 ), ..., e n (y 1 )} is an orthonormal basis, we see that

g jk | γ = δ jk | γ .
This implies that ∂ 1 g jk | γ = 0. Now note that:

∂ α g ij | γ = ∇ g ∂α ∂ i , ∂ j g | γ + ∂ i , ∇ g ∂α ∂ j g | γ ,
where ∇ g denotes the Levi-Civita connection on ( M , g). Using the symmetry for Levi-Civita connection we see that

∇ g ∂α ∂ 1 | γ = ∇ g ∂ 1 ∂ α | γ = ∇ g γ(y 1
) e α (y 1 )| γ = 0, which together with ∇ g γ γ = 0 implies that ∂ j g 11 | γ = 0 for all j ∈ {1, . . . , n}. This implies that Γ 1 1j | γ = 0 for all j, where Γ i jk denotes the Christoffel symbol for ( M , g). Pick an arbitrary unit vector (v 2 , . . . , v n ) ∈ R n-1 and for each y 1 ∈ R consider the geodesic in M parametrized as h(r) = exp γ(y 1 ) (r n α=2 v α e α (y 1 )) with the corresponding Fermi coordinates (y 1 , rv 2 , . . . , rv n ). Note that ḣα (0) = v α , and h satisfies 0 = ḧk αβ (r) + Γ k αβ (r) ḣα (r) ḣβ (r) = Γ k αβ (r) ḣα (r) ḣβ (r), since ḧk αβ (r) = 0 in this coordinate system. Since v α is arbitrary, we deduce that Γ k αβ | γ = 0 for all 1 k n and all 2 α, β n. To conclude that all the Christoffel symbols vanish on the null geodesic, we still need to show Γ α 1β | γ = 0. Using the definition of the Christoffel symbol we see that it suffices to show that

∂ α g 1β | γ = 0. But, ∂ α g 1β | γ = ∂ 1 , ∇ g ∂ α ∂ β g | γ = Γ 1 αβ | γ = 0.
Thus in y coordinates ḡ| β = η jk where η jk denotes the Minkowski metric on R 1+n and ∂ i ḡjk | β = 0. Since the map y → z is affine, it is easy to verify that ḡ satisfies the claimed properties.

It should be clear now that the constants a, a 0 , s -, s + , b 0 , b defined in (3.4) merely denote the s-coordinates of the points β(τ --ε), β(τ --ε 2 ), β(τ -), β(τ + ), β(τ + + ε 2 ), and β(τ + + ε) respectively.

22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE 3.2. Eikonal and Transport equations. Throughout this subsection we will assume that β is a null geodesic in D that is extended to M as described above with local coordinates (z 0 , . . . , z n ) and that A, q satisfy (1.1). Let us consider the differential operator:

L A,q := -∆ ḡ + A∇ ḡ + q. We start the construction of the approximate Gaussian beam by defining the set

V = {(z 0 , z ) ∈ M | z 0 ∈ [a 0 , b 0 ], |z | < δ}, (3.6)
with 0 < δ < δ (see Lemma 3.1) sufficiently small such that the set V does not intersect the sets {0} × M and {T } × M (this can always be fulfilled as β ∈ D). We make a WKB ansatz of the form

u = e iρϕ v, (3.7)
such that u is an approximate solution to L A,q u = 0 as a formal power series in ρ > 1.

Here, ϕ ∈ C 3 (V) and v ∈ C 2 (V). We make the following ansatz for ϕ, v respectively:

ϕ = 2 k=0 ϕ k (s, z ) and v(s, z ) = ρ n 4 v 0 (s)χ( |z | δ ), (3.8)
where for each k = 0, 1, 2, ϕ k is a homogeneous polynomial of degree k with respect to the variables z i with i ∈ {1, ..., n}. The smooth function χ : R → [0, ∞] satisfies χ(t) = 1 for |t| 1 4 , and χ(t) = 0 for |t| 1 2 . We also define the set where Note that: ∆ ḡ(e iρϕ v) = e iρϕ (-ρ 2 dϕ, dϕ ḡv + iρ(2 dϕ, dv ḡ + (∆ ḡϕ)v) + ∆ ḡv).

Let us define:

Sϕ := dϕ, dϕ ḡ, and T A v := 2 dϕ, dv ḡ + (-A∇ ḡϕ + ∆ ḡϕ)v.

Then:

L A,q (e iρϕ v) = e iρϕ (ρ 2 (Sϕ)v -iρT A v + L A,q v). (3.9)
We require that Sϕ = n k,l=0

ḡkl ∂ k ϕ ∂ l ϕ vanishes up to second order on the null geodesic β with respect to the transversal directions (the case Sϕ ≡ 0 is known as an eikonal equation). Put differently, in terms of the Fermi coordinates we require that:

∂ α 1 ∂z 1 α 1 . . . ∂ αn ∂z n αn (Sϕ)(s, 0, . . . , 0) = 0 for s ∈ (a 0 , b 0 ), (3.10)
for all m = 0, 1, 2 and all choices of α 1 , . . . , α n 0 such that n j=1 α j = m. We will also require that the following transport equation holds along the null geodesic β:

(T A v 0 )(s, 0, . . . , 0) = 0 for s ∈ (a 0 , b 0 ). Recalling that ḡ| β = 2dz 0 dz 1 + (dz 2 ) 2 + . . . + (dz n ) 2 , this reduces to

2∂ 0 ϕ ∂ 1 ϕ + n k=2 (∂ k ϕ) 2 = 0. (3.12)
Similarly, for m = 1, we obtain (recall that for all i, j, k we have

∂ i g jk | β = 0) ( n k,l=0 ḡkl ∂ 2 kα ϕ ∂ l ϕ)| β = 0, for all 1 α n. (3.13)
Recalling the definition of the phase function ϕ from equation (3.8), it is clear that equations (3.12) and (3.13) will be satisfied if we set ϕ 0 = 0 and ϕ 1 = r. Next we consider the case m = 2 in equation (3.10) and write ϕ 2 (s, z ) := 1 i,j n H ij (s)z i z j where H ij = H ji is a symmetric matrix. Let us impose the auxiliary condition that

H(s) > 0 for s ∈ (a 0 , b 0 ). (3.14)
This assumption will lead to a Gaussian decay away from the null geodesic β, but we will also provide a geometric motivation behind this assumption in the next section. We require:

∂ 2 ∂z i ∂z j ( n k,l=0 ḡkl ∂ϕ ∂z k ∂ϕ ∂z l )| β = 0 1 ∀i, j n.
This is equivalent to:

(∂ 2 ij ḡkl ∂ k ϕ ∂ l ϕ + 2ḡ kl ∂ 3 kij ϕ ∂ l ϕ + 2ḡ kl ∂ 2 ki ϕ ∂ 2 lj ϕ + 4∂ i ḡkl ∂ 2 jk ϕ ∂ l ϕ)| β = 0. which again simplifies to (∂ 2 ij ḡ11 + 2ḡ 10 ∂ 3 0ij ϕ + 2 n k=2 ∂ 2 ki ϕ ∂ 2 kj ϕ)| β = 0.
We therefore obtain the following Riccati type ODE:

d ds H + HCH + D = 0, s ∈ (a 0 , b 0 ) H(s -) = H 0 with H 0 > 0, (3.15)
where C is the matrix defined through Lemma 3.2. The Riccati equation (3.15) has a unique solution. The solution H is symmetric and (H(s)) > 0 for all s ∈ (a 0 , b 0 ). We have H(s) = Z(s)Y (s) -1 where Z(t) and Y (t) solve the following system of first order linear ODEs:

   C 11 = 0 C ii = 2 2 i n C ij = 0 otherwise
d ds Y = CZ, Y (s -) = I, d ds Z = -DY, Z(s -) = H 0 .
In addition, Y (s) is non-degenerate for all s ∈ [a 0 , b 0 ].

Lemma 3.3. The following identity is satisfied:

det( (H(s)) • | det(Y (s))| 2 = det( (H 0 )).
Let us make some remarks about the regularity of the solutions Y (s), H(s). Since C is a constant matrix, the matrix Y (s) also satisfies

(∆ ḡϕ)| β = N i,j=0 ḡij ∂ 2 ij ϕ| β = n j=2 ∂ 2 jj ϕ| β = Tr(CH).
Thus equation (3.11) simplifies to:

2∂ s v 0 + (Tr(CH) -A(s, 0) β)v 0 = 0 s ∈ [a 0 , b 0 ], (3.19)
where A β := A∂ s = A, dz 1 ḡ. We proceed to prove that

v 0 (s) = det(Y (s)) -1 2 e 1 2 ( s s - A(τ,0) β dτ ) s ∈ [a 0 , b 0 ] (3.20)
Tr(C(s)H(s)) = Tr(C(s)Z(s)Y (s) -1 ) = Tr( dY ds (s)Y (s) -1 ) = d ds log(det(Y (s))),
where we have used the fact that dY ds (s) = C(s)Z(s). Clearly v 0 ∈ C 2 ([a 0 , b 0 ]), which together with the definition of v implies that v ∈ C 2 (V). This concludes the construction of the amplitude function and also the construction of u defined by (3.7) which we refer to as an approximate Gaussian beam.

3.5. Geometrical interpretations. We will briefly discuss some geometrical aspects of the approximate Gaussian beam construction. In particular, we will discuss explicitly, how conjugate points on (M, ḡ) manifest themselves in the vector valued function Y (s) constructed above. First we note the following lemma. Lemma 3.4. Let β be a null geodesic as above and let us consider the Fermi coordinates near β. We have the following identity:

∂ 2 ḡ11 ∂z i ∂z j | β = -2R 0i0j | β ,
for all indices i, j ∈ {1, . . . , n}, where R λµνκ denotes the curvature tensor for the Lorentzian manifold ( M, ḡ).

Proof. First note that by Bianchi identities the expression is clearly symmetric with respect to indices i, j. We let ∇ ḡ and Γi jk denote the Levi-Civita connection and the Christoffel symbol for (M, ḡ) respectively. Recall that ḡij denotes the inverse of the matrix ḡij and since ḡ| β = 2dz 0 dz 1 + (dz 2 ) 2 + . . . + (dz n ) 2 it follows that

∂ 2 ḡ11 ∂z i ∂z j | β = - ∂ 2 ḡ00 ∂z i ∂z j | β .
Since [∂ i , ∂ j ] = 0, we have

R 0i0j | β = ∇ ḡ ∂ 0 ∇ ḡ ∂ i ∂ 0 , ∂ j ḡ| β -∇ ḡ ∂ i ∇ ḡ ∂ 0 ∂ 0 , ∂ j ḡ| β .
Using the definition of the Levi-Civita connection we have

∇ ḡ ∂ 0 ∇ ḡ ∂ i ∂ 0 = n k=0 (∂ 0 Γk i0 )∂ k + n k,m=0 Γk i0 Γm 0k ∂ m , (3.21)
and Recall that ∂ i g jk | β = 0. This implies that Γi jk | β = 0 but since ∂ 0 also denotes the tangent vector to β, we observe additionally that ∂ 0 Γ i jk | β = 0 which implies that (3.21) vanishes along β. Hence we obtain

∇ ḡ ∂ i ∇ ḡ ∂ 0 ∂ 0 = n k=0 (∂ i Γk 00 )∂ k + n k,m=0
R 0i0j | β = - n k=0 ∂ i Γk 00 ∂ k , ∂ j ḡ| β = -( n k=0 ḡjk ∂ i Γk 00 )| β = 1 2 ( n k,m=0 ḡjk ḡkm ∂ 2 im ḡ00 )| β = 1 2 ( n m=0 δ jm ∂ 2 im ḡ00 )| β = 1 2 ∂ 2 ij ḡ00 | β = - 1 2 ∂ 2 ij ḡ11 | β .
We can next use the above lemma in conjunction with the product structure of the Lorentzian manifold M to derive the following corollary.

Corollary 3.5. For any null geodesic β as above, we have the following

∂ 2 ḡ11 ∂z 1 ∂z i | β = 0,
for all indices i ∈ {1, . . . , n}.

This corollary can be used to simplify the Riccati equation (3.15) further. Indeed, Corollary 3.5 implies that D 1i = 0 for i ∈ {1, . . . , n} and since C 1j = C j1 = 0 for all j ∈ {1, . . . , n}, we can simply take H 11 = s + c 0 for any constant c 0 with (c 0 ) > 0, H 1j = H j1 = 0 for all j > 1 and take H i+1,j+1 := Hi,j for all i, j ∈ {1, . . . , n -1} where H is a symmetric (n -1) × (n -1) matrix satisfying

d ds H + 2 H2 + D = 0, H(s -) = H0 , (3.23)
with Dij = D i+1,j+1 for all i, j ∈ {1, . . . , n -1}. This observation also simplifies the construction of the matrix Y . Indeed, we can take Y 11 = c 0 and Y 1j = Y j1 = 0 for all j ∈ {1, . . . , n}.

Note that a key ingredient in the construction of Gaussian beams is the requirement that the matrix valued function Y (s) is non-singular. This is indeed guaranteed in the above construction as a consequence of choosing (H 0 ) > 0. We will briefly discuss what happens when one pursues real valued solutions to this linear system. Recall that Y (s) satisfies equation (3.17). This of course implies that the columns of the matrix Y should also satisfy the same ODE. Let V be one of the columns of Y with representation V = n j=1 V j (s) ∂ ∂z j . Using the definition of the matrix D and Lemma 3.4, we deduce that using the Bianchi identities we obtain (recall that on the null geodeisc, ∂ 0 = β):

d 2 ds 2 V i = n j=1 R 0i0j V j .
D 2 ds 2 V + R(V, β) β = 0
. This is the well-known Jacobi equation along β. We therefore see that the columns of Y are variation fields of some variation of β through null geodesics. In particular based on this geometric characterization of Y , one can deduce that if there exists a point β(s) on the interval [a 0 , b 0 ] that is conjugate to β(a 0 ), then any real valued solution Y (s) to (3.17) will always become singular at that point (see for example [START_REF] Carmo | Riemannian Geometry[END_REF]Section 5.5]). Therefore a global geometric optic construction with a real valued phase function can not be achieved in the presence of conjugate points on (M, ḡ).

3.6. Construction of the remainder terms. With the WKB construction complete, we now return to the task of constructing solutions u 1 , u 2 to (3.1), concentrating on a null geodesic β ∈ D. In particular, we will construct the remainder terms in equations (3.2)-(3.3). We consider the differential operators L A 1 ,q 1 and L * A 2 ,q 2 (formal adjoint of L A 2 ,q 2 with respect to the real L 2 inner product). We can use the previous discussion to obtain two families of approximate solutions given by e iρϕ v 1 and e -iρϕ v2 to these differential operators. Indeed, let

F 1,ρ = -L A 1 ,q 1 (e iρϕ v 1 ), F 2,ρ = -L * A 2 ,q 2 (e -iρ φ v2
). Applying equation (3.9), we obtain

F 1,ρ = -e iρϕ [ρ 2 (Sϕ)v 1 -iρT A 1 v 1 + L A 1 ,q 1 v 1 ], F 2,ρ = -e -iρ φ[ρ 2 (Sϕ)v 2 + iρT -Ā2 v 2 + L * A 2 ,q 2 v2 ].
(3.24)

The phase function ϕ ∈ C 3 (V) is chosen exactly as in Section 3.3 and adapting equation (3.11) to this case, we make the following ansatz for v 1 , v 2 :

v i = ρ n 4 v i,0 (s)χ( |z | δ ) for i = 1, 2,
such the functions v i,0 (s) satisfy the following transport equations:

2∂ s v 1,0 + (Tr (CH) -A 1 β)v 1,0 = 0, s ∈ [a 0 , b 0 ] 2∂ s v 2,0 + (Tr (CH) + Ā2 β)v 2,0 = 0, s ∈ [a 0 , b 0 ]. (3.25)
Using (3.20) we have:

v 1,0 (s) = det(Y (s)) -1 2 e 1 2 ( s s - (A 1 β)(τ,0) dτ ) , v 2,0 (s) = det(Y (s)) -1 2 e -1 2 ( s s - ( Ā2 β)(τ,0) dτ ) . (3.26)
Note that F 1,ρ , F 2,ρ are compactly supported in a small tubular region around the null geodesic where the Fermi coordinates are well defined. Also recall from the previous discussions that v 1 , v 2 ∈ C 2 (V) in the Fermi coordinates. Next we define the expression R j,ρ , j = 1, 2, as the solution of the following IBVP   

L A 1 ,q 1 R 1,ρ = F 1,ρ , (t, x) ∈ M, R 1,ρ (0, x) = 0, ∂ t R 1,ρ (0, x) = 0, x ∈ M R 1,ρ (t, x) = 0, (t, x) ∈ (0, T ) × ∂M, (3.27)
   L * A 2 ,q 2 R 2,ρ = F 2,ρ , (t, x) ∈ M, R 2,ρ (T, x) = 0, ∂ t R 2,ρ (T, x) = 0, x ∈ M R 2,ρ (t, x) = 0, (t, x) ∈ (0, T ) × ∂M. (3.28)
The energy estimate (2.5) in Section 2 implies that equations (3.27) and (3.28) admit unique solutions

R j,ρ ∈ C([0, T ]; H 1 0 (M )) ∩ C 1 ([0, T ]; L 2 (M )) j = 1, 2, with the estimates: R j,ρ H 1 (M) C F j,ρ L 2 (M) . (3.29) We claim that R j,ρ , j = 1, 2, satisfy the following decay property lim ρ→+∞ ( R j,ρ L 2 (M) + ρ -1 R j,ρ H 1 (M) ) = 0, (3.30)
and showing this completes the construction of the solutions u 1 , u 2 of (3.1). Note that for j = 1, 2, using (3.10), (3.11) we have the following bounds:

v j C 2 (M) Cρ n 4 , |T A j v j | Cρ n 4 |z |χ( |z | δ ), |Sϕ| C|z | 2 N (|z |) for z ∈ V, (3.31)
where C > 0 only depending on the geometry and

A 1 C 1 , A 2 C 1 .
Here, N denotes a continuous function depending on the geometry ( M, ḡ) and such that N (0) = 0. Note that equation (3.14) implies that:

|e iρϕ | = |e -iρ φ| e -Dρ|z | 2 , (3.32)
with D > 0 independent of ρ and only depending on the geometry. In particular these estimates imply that: Combining these bounds with (3.24), we find

ρ n 4 e iρϕ 2 L 2 (V) V ρ n 2 e -Dρ|z | 2 χ 2 ( |z | δ )dz = O(1), ρ n 4 (Sϕ)e iρϕ 2 L 2 (V) V |z | 4 ρ n 2 N 2 (|z |)e -Dρ|z | 2 χ 2 ( |z | δ )dz = o(ρ -2 ), (T A j v j )e iρϕ 2 L 2 (V) V ρ n 2 |z | 2 e -Dρ|z | 2 χ 2 ( |z | δ )dz = O(ρ -1 ), (3.33
F j,ρ L 2 (M) = o(ρ), j = 1, 2, (3.34)
and using the estimate (3.29), we deduce that

lim ρ→+∞ ρ -1 R j,ρ H 1 (M) C lim ρ→+∞ ρ -1 F j,ρ L 2 (M) = 0, j = 1, 2.
Therefore, in order to prove the bound (3.30), it only remains to prove that

lim ρ→+∞ R j,ρ L 2 (M) = 0, j = 1, 2. (3.35)
Let us begin by stating the following two lemmas that will simplify the proof of the estimate (3.35).

Lemma 3.6. Let ϕ, V be as above. Then, ∂ t ϕ does not vanish in V.

Proof. Since we are considering the neighborhood V we may use the Fermi coordinate system (s, z ). Recall that in this coordinate system ϕ(z

) = z 1 + H ij (z 0 )z i z j = r + H ij (s)z i z j
where i, j ∈ {1, . . . , n}. Therefore:

√ 2∂ t ϕ = ∂ s ϕ -∂ r ϕ = -1 + 1 i,j n Ḣij (s)z i z j + 2 n i=1 H i1 (s)z i , which implies that for |z | = |z 1 | 2 + . . . + |z n | 2 < δ sufficiently small we have |∂ t ϕ| > 1 2 . (3.36)
Lemma 3.7. Let ϕ be as above and suppose that f ∈ C(V) with supp f ⊂ V. The following estimate holds:

lim ρ→∞ ρ n 4 t 0 f (τ, •)e iρϕ(τ,•) dτ L 2 (M) = 0. (3.37)
Proof. Note that since f is supported near the null geodesic β, we may use the Fermi coordinate system (s, z ). Define

ζ ρ : R n+1 → R through ζ ρ (x) = κ -1 ρ n+1 4 χ(ρ 1 4
|x|), with χ defined as in (3.8) and κ = χ L 1 (R n ) . We define the smooth functions:

f ρ = f * ζ ρ (s, z ) ∈ V.
It is clear that f ρ is supported near the null geodesic β and

lim ρ→∞ f ρ -f C(M) = 0, f ρ W k,∞ (M) C k ρ k 4 ∀k ∈ N. ( 3
ρ n 4 t 0 f (τ, •)e iρϕ(τ,•) dτ L 2 (M) I 1 + I 2 ,
where

I 1 = ρ n 4 t 0 f ρ (τ, •)e iρϕ(τ,•) dτ L 2 (M)
, and

I 2 = ρ n 4 t 0 (f -f ρ )(τ, •)e iρϕ(τ,•) dτ L 2 (M) .
For I 1 , we integrate by part in time using e iρϕ = -i ρ∂tϕ ∂ t e iρϕ and use Lemma 3.6 together with the fact that the set V does not intersect the sets {0}×M and {T }×M to obtain

I 1 ρ n 4 ρ∂ t ϕ f ρ e iρϕ L 2 (M) + ρ n 4 -1 t 0 ∂ τ ( f ρ ∂ τ ϕ )e iρϕ dτ L 2 (M) ρ -3 4 ,
where we are using the first bound in (3.33) together with equation (3.38) to bound the first term by ρ -1 and the second term by ρ -3 4 . For the term I 2 , we can use the Cauchy Schwarz inequality along with inequalities (3.33) and equation (3.38) to obtain

I 2 T ρ n 4 (f -f ρ )e iρϕ L 2 (M) T f ρ -f C(M) • ρ n 4 e iρϕ L 2 (V) → 0 as ρ → ∞.
Proof of Estimate (3.35). The result for R 1,ρ and R 2,ρ being similar, we will only give a proof for R 1,ρ . Using the energy estimate (2.6) in Section 2, together with the fact that R 1,ρ solves equation (3.28), it suffices to prove that lim ρ→∞ F * ,ρ L 2 (M) = 0, where F * ,ρ (t, x) := t 0 F 1,ρ (τ, x) dτ . We have:

F * ,ρ (t, x) = I 1 (t, x) + I 2 (t, x) + I 3 (t, x), (3.39)
where

I 1 (t, x) = - t 0 e iρϕ ρ 2 (Sϕ)(τ, x)v 1 (τ, x) dτ I 2 (t, x) = t 0 e iρϕ [iρT A 1 v 1 (τ, x)] dτ I 3 (t, x) = - t 0 e iρϕ [(L A 1 ,q 1 )v 1 (τ, x)] dτ.
We proceed to bound each of the above integrals using the Fermi coordinates (s, z ) around β. This can be done since each of the above integrands is supported in a small tubular neighborhood of β. For the first integral, I 1 , we apply integration by 22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE parts to obtain

I 1 (t, x) = t 0 iρ 2 ρ∂ τ ϕ (∂ τ e iρϕ )(Sϕ)(τ, x)v 1 (τ, x) dτ = iρe iρϕ ∂ t ϕ (Sϕ)(t, x)v 1 (t, x) + t 0 e iρϕ -iρ∂ τ [( Sϕ ∂ τ ϕ )v 1 ] (τ, x) dτ = iρe iρϕ ∂ t ϕ (Sϕ)v 1 + iρ t 0 e iρϕ ∂ 2 τ ϕ (∂ τ ϕ) 2 (Sϕ)v 1 dt -iρ t 0 e iρϕ ∂ τ [(Sϕ)v 1 ] ∂ τ ϕ dt := S 1 (t, x) + S 2 (t, x) + S 3 (t, x).
(3.40) Using (3.33) we deduce that S 1 L 2 (M) = o(1). For S 2 (t, x), we first use the Cauchy-Schwarz inequality to observe that

S 2 L 2 (M) ρT e iρϕ ∂ 2 t ϕ (∂ t ϕ) 2 (Sϕ)v 1 L 2 (V) .
Recalling that ϕ ∈ C 3 (V) and applying the bound (3.33) again, we deduce that S 2 L 2 (M) = o(1). For S 3 , we integrate by parts again to obtain

S 3 = -e iρϕ ∂ t [(Sϕ)v 1 ] (∂ t ϕ) 2 + t 0 e iρϕ ∂ τ [ ∂ τ ((Sϕ)v 1 ) (∂ τ ϕ)
(M) = o(ρ - 1
2 ) and S 5 L 2 (M) = o(1). Let us now consider the function I 2 (t, x). We have:

I 2 (t, x) = t 0 -i ρ∂ τ ϕ ∂ τ (e iρϕ ) [iρT A 1 v 1 (τ, x)] dτ = e iρϕ ∂ t ϕ (T A 1 v 1 )(t, x) - t 0 e iρϕ ∂ τ ( T A 1 v 1 ∂ τ ϕ ) dτ := S 6 (t, x) + S 7 (t, x).
(3.41) Using (3.33) again, we deduce that S 6 L 2 (M) ρ -1 2 . For S 7 , noting that ∂ t (

T A 1 v 1 ∂tϕ ) ∈ C c (V)
, we can use Lemma 3.7 to show that lim ρ→∞ S 7 L 2 (M) = 0. Finally, for the function I 3 (t, x), since ϕ, v 1 ∈ C 2 (V), we have L A 1 ,q 1 v 1 ∈ C(V) and that this function is supported near the null geodesic β. Using this observation together with Lemma 3.7, we conclude that lim ρ→∞ I 3 L 2 (M) = 0. This completes the proof.

ρ n 2 e -2ρ (ϕ) (det |Y (s)|) -1 χ 2 ( |z | δ ) dz = C + O(|ρ| -∞ ),
where C = 0 is a constant. Combining Fubini's theorem and the above equation we have

0 = b 0 a 0 Ġ(s) ds = s + s - Ġ(s) ds = G(s + ) -G(s -).
Since G(s -) = 1 by definition, we conclude that G(s + ) = 1, which implies that for any null geodesic β ∈ D, we have

exp( 1 2 β A) = 1.
Since supp A ⊂ E, we can conclude that:

1 4πi β A ∈ Z, (4.3)
for any maximal null geodesic in R × M . Now consider equation (4.3) and write β(t) = ( t + t, γ(t)), where t ∈ I and t ∈ R. We consider the family of null geodesics β s (t) = (s + t, γ(t)) and note that equation (4.3) holds for all s ∈ R. Since A is of compact support, we must have that βs A vanishes for |s| large. Together with the continuity of A, we can conclude that it must vanish for all s. Therefore, for any maximal null geodesic in R × M we have

β A = 0.
Finally, under the hypothesis of Theorem 1.3, together with Proposition 1.4, we conclude that the first equation in (1.4) holds.

4.2.

Recovery of the light ray transform of q. The previous discussion yields that A 1 -A 2 = dψ for some ψ ∈ C 2 (M) with ψ| ∂M = 0. We define ψ = 1 2 ψ, and let The gauge invariance of the DN map implies that

Ã2 (t, x) = A 2 (t, x) + 2d ψ = A 1 (t, x) ∀(t, x) ∈ M, q2 (t, x) = q 2 (t, x) + ∆ ḡ ψ -A 2 ∇ ḡ ψ -∇ ḡ ψ, ∇ ḡ ψ ḡ ∀(t, x) ∈ M.
Λ A 1 ,q 1 = Λ Ã2 ,q 2 = Λ A 1 ,q 2 . (4.5)
Furthermore, we also know that supp (q 2 -q 1 ) ⊂ E. We now proceed as in the previous section. We start by fixing q = q 1 -q2 . As before, we take β to be a maximal null geodesic in D ⊂ M and extend it to M. For j = 1, 2, we define

u j ∈ C([0, T ]; H 1 (M )) ∩ C 1 ([0, T ]; L 2 (M )),
as the Gaussian beam solutions of (3.1) (with A 2 replaced with Ã2 and q 2 replaced with q2 ), taking the form (3.2)-(3.3). We fix u 3 the solution of (1.3) with A = Ã2 , q = q2 and f given by

f (t, x) := u 1 (t, x) = v 1 (t, x)e iρϕ , (t, x) ∈ (0, T ) × ∂M. Then the function u = u 3 -u 1 solves    -∆ ḡu + Ã2 (t, x)∇ ḡu + q2 (t, x)u = qu 1 , in M, u = 0, on (0, T ) × ∂M, u(0, •) = 0, ∂ t u(0, •) = 0 in M. (4.6)
Moreover, (4.5) implies

∂ ν u = ∂ ν u 3 -∂ ν u 1 = Λ Ã2 ,q 2 f -Λ A 1 ,q 1 f = 0. Multiplying ( 4
-2ρ ϕ v 1 v2 dV ḡ = 0, (4.8)
Here, we note that since Ã2 = A 1 , we have

v 1 v2 = ρ n 2 v 1,0 v2,0 χ 2 ( |z | δ ) = ρ n 2 det(|Y (s)|) -1 χ 2 ( |z | δ ).
Thus, taking the limit ρ → ∞ and using Lemma 3.3 we deduce that t∈I q(β(t)) dt = 0. (4.9)

This equation only holds for maximal null geodesics in D, but since supp q ⊂ E, we can conclude that this equation holds for all maximal null geodesics β in R × M . Together with Proposition 1.4, we conclude that q = 0, thus completing the proof of Theorem 1.3.

22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE

Inversion of the light ray transforms

This section is concerned with the proof of Proposition 1.4. We start with the inversion of the light ray transform of a scalar function q satisfying supp q ⊂ E.

Proof of statement (i) in Proposition 1.4. We know that for any maximal geodesic

β in R × M , t∈I q(β(t)) dt = 0.
Using the identification of maximal null geodesics β = ( t + t, γ(t)) with maximal geodesics in M , we conclude that

Lq( t, γ) = I q( t + t, γ(t)) dt = 0, ( 5.1)
for all t ∈ R and all unit speed maximal geodesics γ in M . Taking the Fourier transform of L( t, γ) with respect to the variable t ∈ R and using the fact that q is compactly supported, we deduce that

0 = Lq(τ, γ) = I R e -iτ tq(t + t, γ(t)) d t dt = I e iτ t q(τ, γ(t)) dt,
where q(τ, •) denotes the Fourier transform of q with respect to the time variable.

Evaluating at τ = 0, we get I q(0, γ(t)) dt = 0, for all geodesics γ in M . Using Hypothesis 1.2, we deduce that q(0, •) = 0. Next, we evaluate all the derivatives of L(τ, γ) at τ = 0:

0 = [∂ k τ Lq(τ, γ)]| τ =0 = [∂ k τ I e iτ t q(τ, γ(t)) dt]| τ =0 = k j=0 k j I (it) k-j ∂ j τ q(0, γ(t)) dt.
(5.2) We now use induction on j to show that ∂ j τ q(0, •) = 0 for all j ∈ N. This is true for j = 0 as shown above. Suppose the induction hypothesis holds for all j < k. Then using (5.2) implies that

I ∂ k τ q(0, γ(t)) dt = 0.
Together with Hypothesis 1.2, we conclude that ∂ k τ q(0, •) = 0. This completes the proof by induction. Now, since q(t, •) is a compactly supported function in t, we know that q(τ, •) is real analytic with respect to τ and since all the derivatives vanish at τ = 0, we conclude that q vanishes identically, which implies that q ≡ 0.

22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE It remains to prove the inversion of the light ray transform of a one-form A satisfying supp A ⊂ E up to the natural gauge. let us begin with some remarks and lemmas. We write A = b dt + ω. Similar to the proof above, we have

0 = LA(τ, γ) = R e -iτ t I [b( t + t, γ(t)) + ω( t + t, γ(t)) γ(t)] dt d t,
for all maximal geodesics γ ∈ M . Since A is compactly supported, we can interchange the order of integration to obtain

0 = LA(τ, γ) = I ( b(τ, γ(t)) + ω(τ, γ(t)) γ(t)) dt,
where b and ω denote the Fourier transform in time of the compactly supported function b and the one-form ω respectively. Evaluating at τ = 0, we obtain

I γ ( b(0, γ) + ω(0, γ)) = 0.
Using Hypothesis 1.2 together with smoothness properties of the Helmholtz decomposition discussed in Section 1.3, we deduce that b(0, x) = 0, ω(0, x) = dψ 0 (x), (

for some ψ 0 ∈ C 1 (M ) vanishing on ∂M . Note that since ω is C 1 smooth, and since dψ 0 (x) = ω(0, x) ∈ C 1 (M ; T * M ), it follows that

ψ 0 ∈ C 2 (M ).
We have the following lemma.

Lemma 5.1. Let ψ -1 := 0 and ψ 0 be as above. There exists a sequence of functions {ψ k } ∞ k=1 ⊂ C 2 (M ) all vanishing on ∂M and such that for every k 0,

∂ k τ ω(0, x) = dψ k (x), and ∂ k τ b(0, x) = ikψ k-1 (x). Proof.
We proceed by induction. The claim is true for k = 0 as discussed above. Suppose it holds for all m < k. For m = k, we have: 22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE Now using the fact that I t k-j dψ j γ(t) dt = -I (k -j)t k-j-1 ψ j dt, we can simplify the right hand side of the above equation to Note that ψ| (0,T )×∂M = 0 as ψ k | ∂M = 0 for all k 0. We will now show that ψ ∈ C 2 (M). Indeed, it is clear from equation (5.4) that ψ ∈ C 1 (M). But then since dψ = A ∈ C 1 (M; T * M) we may conclude that ψ ∈ C 2 (M). 5.1. Remarks. With the proof of Theorem 1.3 complete, let us state a few remarks. We start with the auxiliary condition that the coefficients A, q should be known on the set M \ E. This is merely an artifact of the light ray inversion method, as the Gaussian beam construction shows that the light ray transforms of A and q can be obtained on the set D. As the ratio D g (M )/T grows, the set E grows and the coefficients are therefore recovered in a large set that is closer in size to the optimal set D. The assumption that D g (M ) < ∞ in particular implies that the manifold (M, g) should be non-trapping. To illustrate this remark, consider M = [0, 1] × S 1 with S 1 denoting the unit circle and note that in this case D g (M ) = ∞.

Regarding the smoothness of the metric g, note that we only require the metric in Fermi coordinates to be C 4 smooth. This however, can only be guaranteed if the 22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE metric is a priori known to be C 6 smooth, as there could be some loss of regularity in the angular directions of the exponential map of a Riemannian manifold. One can in fact use a mollification of the phase function ϕ to improve the result to g ∈ C 4 (M ; Sym 2 M ). The key here is that there is no loss of regularity in the direction tangent to the null geodesic in Fermi coordinates. We believe that the result could be extended to g ∈ C 2 (M ; Sym 2 M ), but this will require mollification of the metric, g ρ . Improvements beyond the C 2 -smoothness for the metric g could be much harder.

3 . 3 .

 33 Construction of the phase. We begin by solving equation(3.10). For m = 0, we obtain the equation ( β = 0 1 ∀i, j n.

1 -

 1 Submitted to Analysis & PDE and D = (D ij ) n×n := 1 4 (∂ 2 ij ḡ11 | β ) n×n for all i, j ∈ {1, . . . , n}. Note that since g is C 4 smooth in the Fermi coordinates, we have D ∈ C 2 ([a 0 , b 0 ]; C n×n ). We now recall two lemmas. For the proofs, we refer the reader to [21, Lemma 8,Section 8] and [21, Lemma 10,Section 8] respectively.

d 2 ds 2 Y

 2 = -CDY, Y (s -) = I, Ẏ (s -) = CH 0 . (3.17) since D ∈ C 2 ([a 0 , b 0 ]; C 2n ), we immediately deduce that Y ∈ C 4 ([a 0 , b 0 ]; C n×n ). Now considering the ODE for the function Z(s), we deduce that Z ∈ C 3 ([a 0 , b 0 ]; C n×n ). Finally, since H(s) = Z(s)Y -1 (s) and since Y (s) is non-singular on [a 0 , b 0 ], we conclude that H ∈ C 3 ([a 0 , b 0 ]; C n×n ) and ϕ ∈ C 3 (V) (3.18) in the Fermi coordinates. 3.4. Construction of the amplitude. Let us now study the transport equation (3.11). First observe that in Fermi coordinates

Rearranging the indices and 22

 22

1 - 2 Ġ

 12 Submitted to Analysis & PDE Using Lemma 3.1 we have || det ḡ|(s, z ) -1| C|z | 2 . Using the fact that A, q = 0 on the set M \ M and the above two obsevations we get lim (s)e -2ρ (ϕ) (det |Y (s)|) -1 χ 2 ( |z | δ ) ds ∧ dz = 0. Lemma 3.3 implies that |z |<δ

. 6 )

 6 by u 2 and integrating by parts we obtain M qu 1 u 2 dV ḡ = 0

  k-j [∂ j τ b(0, γ(t)) + ∂ j τ ω(0, γ(t)) γ(t)] dt = 0.Now using the induction hypothesis together we can simplify this expression to obtainI [∂ k τ b(0, γ(t)) + ∂ k τ ω(0, γ(t)) γ(t)] dt = -k-j [(ij)ψ j-1 + dψ j γ(t)] dt).

- 1 - 6 )

 16 ) k-j (ij)ψ j-1 -(it) k-j-1 i(k -j)ψ j ] dt) ) k-j (ij)ψ j-1 ] dt) + )(it) k-j ψ j-1 ] dt) = ik I ψ k-1 (γ(t)) dt.Hence we haveI [∂ k τ b(0, γ(t)) -ikψ k-1 (γ(t)) + ∂ k τ ω(0, γ(t)) γ(t)] dt = 0.Using Hypothesis 1.2 together with smoothness properties of the Helmholtz decomposition in Section 1.3, we conclude that there existsψ k ∈ C 1 (M ) vanishing on ∂M such that ∂ k τ ω(0, x) = dψ k (x), and ∂ k τ b(0, x) = ikψ k-1 (x). Since dψ k = ∂ k τ ω(0, x) ∈ C 1 (M ; T * M ) and ψ k ∈ C 1 (M ), we conclude that ψ k ∈ C 2 (M ).This completes the induction argument.Proof of statement (ii) in Proposition 1.4. Let us defineψ(t, x) := t 0 b(s, x) ds = t -∞ b(s, x) ds,(5.4)where we are using the fact that supp A ⊂ M. Note that since A has compact support, equation (5.3) impliesψ(t, x) = ψ(T, x) = T 0 b(s, x) ds = b(0, x) = 0 for t T . (5.5)Similarly, ψ(t, x) = 0 for t 0. Again, we let ψ(τ, x) denote the Fourier transform of ψ in the time variable. Note that since ψ is compactly supported in time, ψ(τ, x) is analytic with respect to τ . Let us define the coefficients { ψk } ∞ Submitted to Analysis & PDEWe claim that ψk (x) = ψ k (x) holds for all k 0 and all x ∈ M . To see this note that by definition ∂ t ψ = b. Hence iτ ψ(τ ) = b(τ ). Now differentiating this expression k + 1 times and evaluating at τ = 0 we deduce thati ψk = i ψ(k) (0) = i 1 k + 1 b(k+1) (0) = iψ k ,where we used Lemma 5.1 in the last step. Thus, we haveψ(τ, x) = ∞ k=0 ψ k (x) k! τ k .Since ω(t, x) is also compactly supported in t, Lemma 5.1 implies that ω(τ, x) = Note that for every fixed x ∈ M , the following estimate holds:|∂ k τ ω(0, x)| ∂ k τ ω(τ, x) L ∞ (R) T 0 |t| k |ω(t, x)| dt,which implies that the infinite series (5.6) is uniformly convergent with respect to x ∈ M , and therefore we can write ω(τ, x) = k ) = d ψ(τ, x).Hence, ω = dψ, and subsequently A = b dt + ω = (∂ t ψ) dt + dψ = dψ.(5.7)

  ∈ C 6 (M ; Sym 2 M ), we deduce that the pull back of the metric metric ḡ, Φ * ḡ is C 4 smooth in the Fermi coordinate system.

22 Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE g

  2 ] dτ := S 4 + S 5 . Using equation (3.10), it is clear that |∂ t (Sϕ)| |z |N (|z |) and |∂ 2 t (Sϕ)| N (|z |) on the set V. Using this observation, together with bounds analogous to (3.33), we deduce that S 4 L 2

(t + y 1 ) + a 2 , r := z 1 := 1

Jan 2019 11:09:28 PST Version 1 -Submitted to Analysis & PDE

Acknowledgments

A.F was supported by EPSRC grant EP/P01593X/1, J.I. was supported by the Academy of Finland (decision 295853), Y.K. was partially supported by the Agence Nationale de la Recherche grant ANR-17-CE40-0029 and L.O was supported by the EPSRC grants EP/P01593X/1 and EP/R002207/1.

Reduction to the light ray transform

In this section we will obtain the light ray transforms of A and q on null geodesics from the knowledge of the Dirichlet to Neumann map Λ A,q and then use Proposition 1.4 to prove Theorem 1.3. Proposition 1.4 will subsequently be proved in the next section.

4.1.

Recovery of the light-ray transform of A. From now on we fix q = q 1 -q 2 , A = A 1 -A 2 extended by zero to (R × M ) \ M. We take β to be a maximal null geodesic in D ⊂ M and extend it to M. For j = 1, 2, we define

a solution of (3.1) taking the form (3.2)-(3.3) with the properties described in the previous section. We fix u 3 the solution of (1.3) with A = A 2 , q = q 2 and f given by

Moreover, since A 1 = A 2 on (0, T )×∂M by assumption, we have

Recall that the functions v 1 and v 2 are supported in a small tubular neighborhood of the null geodesic β. Thus the integrand in the above equation is supported near β and we can the use Fermi coordinates z = (z 0 , . . . , z n ) = (s, z ) = (s, r, z 2 , . . . , z n ) = (s, r, z ) to compute the limit. We have:

where G(s) = e