Yasmine Mebrouk
email: mebrouk@unice.fr

Franck Mady

Mourad Benabdesselam

Jean-Bernard Duchez

Wilfried Blanc

Experimental evidence of Er 3+ ions reduction in the radiation-induced degradation of erbium-doped silica fibers

Keywords:

The gain of erbium doped fiber amplifiers is damaged by irradiation partly because of creation of color centers responsible of excess absorption at pump and signal wavelengths. Based on the combination of thermally stimulated luminescence and spectrophotometry, this letter demonstrates that that a part of the gain loss should be associated with the reduction of the density of Er 3+ ions by irradiation.

Erbium-doped fibers (EDF) have proved to be very attractive for a large domain of applications such as superluminescent sources, lasers, or optical signals amplifiers emitting around 1.55 µm. In recent years, EDF have been thought to be employed in harsh environments (nuclear reactor, space communications,... ) [1]. However, ionizing radiations damage EDF and lead to the degradation of gain and signal transmission [1,2]. It is admitted that the gain degradation is due to the creation of color centers (CC) that are responsible for a so-called radiation-induced attenuation (RIA), that is excess loss at signal and pump wavelength notably [3][4][5]. If the absorption and emission cross sections of Er 3+ ions are not modified by irradiation [4], these degradations can be attributed to the combination of two effects: (i) pump absorption by CC [3] reduces the pump power available for population inversion therefore the gain is damaged. (ii) Signal absorption by CC. For instance, phosphorus doping introduces absorbing defects near 0.8 eV (1550 nm) that leads to a decrease of the transmitted signal power during irradiation. Aside from the role of CC, the direct involvement of Er 3+ ions in the gain degradation has not been investigated so far, probably because it has been ruled out based on hasty conclusions [4][5][6]. However, a recent work by Likhachev et al. [7] concluded indirectly that Er 3+ ions may contribute to the gain degradation by comparing two germanosilicates (irradiated in identical conditions) with the same amount of GeO 2 but two different concentrations of Er 2 O 3 . These authors found a higher RIA level in the fiber containing more erbium. Based on the coupling between thermally stimulated luminescence (TSL) and spectrophotometry measurements, we demonstrate in a much more straightforward way that the gain degradation is partly due to the reduction of the density of Er 3+ ions.

Samples consisted of preform discs with 0.5 mm core radius and 1 mm thickness and also of corresponding fibers. Preforms were fabricated by the conventional modified chemical vapour deposition (MCVD) process and solution doping technique and drawn in our laboratory. Measurements in the UV-visible were carried out on preform samples to avoid multimodal propagation artifacts in fibers. Near infra-red (NIR) measurements were made on fibers with an Anritsu optical spectrum analyzer. We investigated three silica-based preform compositions: one with erbium only (pEr), two with different concentrations of aluminum and the same amount of erbium (pAl1:Er, pAl2:Er), and one with germanium and erbium (pGe:Er). Density of Er 3+ ions is about 1.5×10 25 m -3 for all samples. pAl1:Er and pAl2:Er contain 0.25 and 1.37 at.% of Al, respectively. pGe:Er contains 0.15 at.% of Ge. Irradiations were carried out with 45 kV x-rays from a copper-anode tube at 1.2 krad min -1 (SiO 2 ). Absorption measurements were performed by a PerkinElmer Lambda 1050 spectrophotometer. Both irradiations and absorption measurements were conducted at room temperature (RT). TSL of irradiated samples was measured by a Harshaw TLD 3500 reader. TSL reveals all the optically active centers and defects whose valence state has been modified by irradiation. It consists of photons emitted upon radiative recombination of "detrapped" carriers with carriers of opposite polarity on reduced/oxidized centers (recombination centers). This standard interpretation framework of TSL, based on trapping and recombination mechanisms, has been shown to be a proper basis to explain the RIA build-up and bleaching [8,9]. After radiation exposure, samples are heated up linearly (2 ˚C s -1 ) from RT up to 600 ˚C and carriers are thermally released from traps of increasing depths. The TSL is plotted as a function of temperature to give "glow curves" that reflect the energy distribution of trapped states generated by irradiation. Glow curves exhibit "peaks". Each trapping level gives rise to a single peak (the deeper the trap, the higher the peak temperature). TSL spectra acquired at each temperature with an optical multi-channel analyzer allow the identification of photoionized luminescent centers.

Normalized TSL glow curves obtained after x-ray irradiation (300 krad(SiO 2 )) of our samples are shown in Fig. 1. They all present three main broad peaks indicating that three groups of traps are populated by irradiation. Peak 1 lays between 50 and 150 ˚C, peak 2 and 3 are in the 200-300 ˚C and 350-450 ˚C range respectively. For the pEr and pAl1:Er samples, peak 2 is not intense enough to be separated clearly from the broad first peak. A similar peak at 250 ˚is found in pGe:Er but in this case it is attributed to the recombination of released electrons on ionized germanium lone pair center (GLPC) [10][11][12]. Fig. 2 represents the spectral analysis of the TSL emissions. Spectra were measured at 100 ˚C but the same luminescence features are observed throughout the temperature range. All samples show green emission bands at 525 and 547 nm corresponding respectively to transitions from 2 H 11/2 and 4 S 3/2 to the ground state level 4 I 15/2 of Er 3+ ions. In the case of the Ge-doped sample (pGe:Er, inset), one can see the well-known intense blue-violet emission of GLPC recombination [10,13]. In addition, one can distinguish a weak Er 3+ emission at 404 and 474 nm in the case of pAl:Er and pEr samples. These features probably exist in the spectrum of pGe:Er but might be hidden by the very intense emission of GLPC. Trap emptying during the TSL readout is followed by recombination of released carriers, which is accompanied by emissions of Fig. 2. Since this TSL mechanism leads to the thermal bleaching of the samples (see below), it must consist of the reverse processes to those operating during irradiation, i.e. ionization leading to generation of free carriers and traps filling responsible for darkening and RIA. The presence of the emission of Er 3+ ions indicates that carriers released from TSL trap recombine with previously ionized centers to form an excited Er 3+ ions that subsequently de-excites. Somehow, this means that erbium ions have been ionized by irradiation and thus that Er 3+ ions have been lost.

To correlate the TSL glow curves with absorption measurements, thermal bleaching of RIA was examined according to the following procedure: after the optical absorption of the pristine preforms core was measured, samples were irradiated with x-rays (300 krad (SiO 2 )). The absorption of the irradiated core was measured and then the sample was heated up to T stop = 50 ˚C at 2 ˚C s -1 (the same heating rate as in TSL readouts to allow comparison). The new core absorption was measured after cooling down at room temperature. Then, the sample was heated up to T stop = 100˚C and the new core absorption was measured, etc ... Al-OHC at 540 nm and NBOHC at 620 nm [15]). At the 378 and 521 nm Er 3+ absorption wavelengths we notice two "holes" corresponding to a decrease of Er 3+ ions absorption by irradiation. The depth of these two "holes" depends on the codoping element. It was found to increase with the Al concentration but is less pronounced in the Ge-doped sample. Since thermal bleaching experiments are done at the same heating rate as TSL readouts, the Gaussian deconvolution of absorption spectra at each T stop (Fig. 3) may allow us to compare the temperatureresolved bleaching of each gaussian component with TSL glow curves and thus to assign related defects to TSL peaks and traps. This task is actually very hard, due to the multiplicity of the deconvolution solutions notably. We analyzed the variation of Er 3+ absorption (α(Er 3+ )) at 378 and 521 nm by extracting the corresponding bands from spectra of Fig. 3 using the best-suited interpolation of each data spectrum; the obtained result is shown in Fig 5-a. Clearly, the 378 and 521 nm bands areas are decreased by irradiation. Then, absorption recovers gradually with increasing T stop . The change of α can be due to a variation either in the absorption cross-section σ abs or in the density of erbium ions N . It has been argued that Er 3+ ions absorption and emission cross-sections are not affected by γ-radiations [4]. Our results also support this finding in case of x-ray irradiation. Indeed, normalized curves of Fig 5-b show that all spectra exhibit the same shape with identical FWHM (0.03 eV at 521 nm and 0.04 eV at 378 nm) before and after irradiation and throughout the thermal recovery protocol. Measurements between 1400 and 1600 nm on Al:Er-doped preforms and fibers (made in our laboratory or commercial), reported in Fig. 6, confirm the decrease of α(Er 3+ ) observed in UV-visible range. The absorption level is reduced after irradiation and normalization of spectra in insets again shows that the band shape does not change after irradiation. Therefore, we conclude that the absorption cross section σ abs does not change and the decrease in α(Er 3+ ) after x-ray irradiation is due to the decrease of the density N of Er 3+ ions. One can underline that the drawing process has no impact on this result. The direct contribution of Er 3+ ions in the gain degradation has been dismissed to date [4,6]. Obviously a significant amount of Er 3+ ions is lost by irradiation. Fig. 7 illustrates the correlation between the TSL curves and the recovery curves of Er 3+ absorption lines for pAl1 (and 2):Er samples. Although not reported here, similar behaviors have been observed for other samples. The ordinate of each point corresponds to the percentage recovery in the 50 ˚C interval. For example: the point at 325 ˚C corresponds to the difference between area under curves at 350 ˚C and 300 ˚C normalized by the RIA at the given wavelength. Recovery plots at 378 and 521 nm have been deduced from Fig. 5. They should be identical.

In practice, they roughly show similar trends. Discrepancies are mainly due to the weak recovery within each 50 ˚C slot, so differences between spectra at two successive T stop are highly inaccurate. The recovery of Er 3+ ions follows the trap emptying depicted by the TSL curve, except for the first peak component below 100 ˚C which is apparently poorly involved in the process (carriers released from the corresponding traps do not strongly take part in the re-formation of Er 3+ ions so additional TSL emissions probably exist out of the wavelength window of our spectral measurements (Fig. 2). A possible explanation of the reduction of Er 3+ density after irradiation is a valence change to divalent ions. The change in erbium ionization state is as likely as in the case of other rare earth ions. Europium and ytterbium are known to exist in divalent state in oxide glasses; reduction of these rare earth ions by irradiation is possible [8,16]. Divalent erbium Er 2+ was observed in CaF 2 [17] and has been obtained by synthetic methods [18] but it has never been reported in an oxide glass. The radiation induced attenuation mechanism in Er-doped silica could be the same as in Yb-doped [8]. According to TSL and absorption measurements, we can propose the following mechanism: under irradiation Er 3+ ions are reduced to the form Er 2+ by trapping an electron, a hole is consequently captured to form trapped hole centers (AlE', Al-OHC, NBOHC [13,14]) that obviously contribute to RIA (Fig 4). Upon thermal or optical excitation, radiation induced trapped holes are released and recombine to reform Er 3+ ions. Based on the variation in Er 3+ ions absorption bands, the density of ions reduced after irradiation (300 krad) is about 20% for pAl1:Er, 32% for pAl2:Er (with more Al) and 12% for pGe:Er sample. In the case of Ge-codoping, the reduction is smaller for the same density of Er 3+ ions. It could be explained as follows: Ge introduces a huge amount of electron traps (Ge (1) and Ge(2)) [13] that are in competition with Er 3+ for trapping electrons and could trap more electrons than Er 3+ . By contrast, Al-doping introduces hole traps and favors electrons trapping on Er 3+ ions. Moreover, Al and Er ions are well-known to be closely located which may favor the recombination process [19]. In small signal conditions, the gain loss resulting from the conversion of erbium ions should be proportional to the conversion rate, say in the 10-35 % range (at 300 Gy) according to our data. For 30 dB amplifier, e.g., the gain should be decreased by 3-10 dB due to the radiation effect on erbium ions (out of RIA). This should be compared to typical RIA levels reported in literature, i.e. about 20 dB for a 2m-long conventional aluminosilicate amplifier at 1550 nm and about 0.4 dB for radiation-hardened fibers, in same conditions [7,20]. Therefore the loss in Er 3+ ions can have a non-negligible impact on the gain degradation, in addition to RIA. This contribution will be more critical in fibers exhibiting low RIA levels due to optimization of the core compositions [7,20]. It could even determine the gain degradation in highly RIA-resistant fibers. Complementary characterizations are of course needed to estimate the consequences of the basic result demonstrated in this letter.

FIG. 1 .

 1 FIG. 1. Normalized TSL intensity as function of temperature for pAl1:Er, pAl2:Er, pEr preform samples. Inset: pGe:Er sample.

2 FIG. 4 .

 24 FIG. 4. Radiation induced absorption (difference between post-and pre-radiation absorption data) for pAl1:Er and pAl2:Er sample and for pEr and pGe:Er in the inset.

FIG. 5 .FIG. 6 .

 56 FIG. 5. Extracted 521 and 378 nm Er 3+ absorption band : (a) before and after irradiation and for all Tstop of thermal recovery. (b): The same normalized with respect to the maximum.

FIG. 7 .

 7 FIG. 7. Thermal recovery of the absorption for pAl1:Er and pAl2:Er preforms at 378 and 521 nm. Normalized TSL curves are also shown for comparison.

  The partial core bleaching within 50 ˚C slots was monitored by considering T stop

		350						500	pGe:Er
		325				4 I 15/2 2 H 11/2	TSL (a.u.)	300 400	Er 3+ GLPC emission emission
	TSL (u.a.)	300		2 H 9/2	2 K 15/2		300 4 I 15/2 4 S 3/2	400	500 nm)	600
		275		4 I 15/2	4 I 13/2			pAl2:Er
						pEr	pAl1:Er
		300	400		500			600	700
						(nm)
	FIG. 2. Normalized TSL spectra after x-rays irradiation (300
	krad) at 100˚C for Al:Er-, Er-and Ge:Er-(inset) doped core
	preforms.
		7	pAl1:Er
		6		2 G
				11/2
		5		378 nm	4 S
						3/2
	(cm -1 )	4				521 nm
		3
									After irradiation
		2
		1	Heating
		0						Before irradiation
		200	300	400	500		600	700	800
						(nm)
									378 and
	521 nm, respectively. Core absorption increases in the
	UV-visible range after irradiation. After heating up to
	600 ˚C we observe the total recovery of RIA and TSL
	readouts do not provide any signal. Therefore the total
	emptying of TSL traps well corresponds to the complete
	RIA bleaching. Absorption by other samples (not shown
	here) exhibits the same darkening/bleaching behavior.
	The difference between the post-and pre-irradiation ab-
	sorption spectra (i.e. RIA) is plotted in Fig 4. RIA plots
	reveal the formation of absorbing centers characterized
	by overlapping absorption bands lying from UV to visi-
	ble (for instance, AlE' absorption band at 300 nm [14],

FIG.

3

. Absorption spectra before and after x-rays irradiation (300 krad) for Al:Er-doped core preforms. Dashed curves : absorption after heating up from RT up to 600 ˚C by steps of 50 ˚C.

= 150, 200, ... up to 600 ˚C. Absorption coefficients α before and after x-ray irradiation of pAl1:Er core sample are plotted in Fig.

3

. Dashed curves correspond to absorption spectra. Er 3+ ground state absorption transitions to 2 G 11/2 and 4 S 3/2 levels are located at

This work, including the PhD thesis of Yasmine ME-BROUK, has been funded by the AirbusGroup foundation under grant # 098-AO11-1108. The authors gratefully acknowledge this support.


Y Grisel

V Mouysset

P A Mazet

J P Raymond

Adaptive refinement and selection process through defect localization for reconstructing an inhomogeneous refraction index

We consider the iterative reconstruction of both the internal geometry and the values of an inhomogeneous acoustic refraction index through a piecewise constant approximation. In this context, we propose two enhancements intended to reduce the number of parameters to reconstruct, while preserving accuracy. This is achieved through the use of geometrical informations obtained from a previously developed defect localization method. The first enhancement consists in a preliminary selection of relevant parameters, while the second one is an adaptive refinement to enhance precision with a low number of parameters. Each of them is numerically illustrated.

Introduction

We are interested in the inverse medium problem consisting in the reconstruction of an inhomogeneous acoustic refraction index from far-field measurements generated through plane waves. This parameter identification problem is non-linear and ill-posed, and we investigate two methods to reduce the number of computed parameters while preserving the reconstruction accuracy. Applications are, for example, non-destructive structure testing or biomedical imaging [START_REF] Scott | A review of non-destructive testing of composite materials[END_REF][START_REF] Song | Through-wall imaging (twi) by radar: 2-d tomographic results and analyses[END_REF][START_REF] Mojabi | Microwave biomedical imaging using the multiplicative regularized gauss-newton inversion[END_REF].

Following the abundant literature, we write the inverse medium problem as a least-squares problem (see [START_REF] Bakushinsky | Iterative methods for approximate solution of inverse problems[END_REF] and references therein). Besides, since we consider discontinuous indices, we look for the index of refraction as a piecewise constant function. In this setting, for its ease of implementation and its efficiency for reasonably sized problems, we consider the Gauss-Newton method, applied to a cost functional involving a Tikhonov regularization [START_REF] Engl | Regularization of inverse problems[END_REF]. However, the Gauss-Newton method treats all parameters in the same way. Yet, during the reconstruction, or even right from the beginning, the values of some parameters can be close to the exact value, while other parameters will need more iterations before reaching a given accuracy. In the absence of some local information, the accurate parameters are then uselessly updated at each iteration. Thus, we explore two uses of geometrical informations, obtained through defect localization, to focus the reconstruction and consequently lighten its numerical cost.

By defect localization, we mean localizing the support of a perturbation with respect to some known reference. However, in this paper, it is the reconstructed index that we use as the known reference, and it is the exact index that we use as an unknown perturbed state. Thus, defect localization can be used to locate errors in the index reconstruction. Besides, it has recently been proved that the location of the defects in a given refraction index could be established from farfield measurements of the unknown state and computed through a fast numerical method [START_REF] Grisel | Determining the shape of defects in non-absorbing inhomogeneous media from far-field measurements[END_REF][START_REF] Bondarenko | The factorization method for inverse acoustic scattering in a layered medium[END_REF]. Also, shape reconstruction has already been used jointly with parameter identification by using the Linear Sampling method [START_REF] Kress | A hybrid method for sound-hard obstacle reconstruction[END_REF][START_REF] Brignone | The use of the linear sampling method for obtaining super-resolution effects in Born approximation[END_REF][START_REF] Brignone | A hybrid approach to 3d microwave imaging by using linear sampling and ACO[END_REF][START_REF] Giorgi | Application of the inhomogeneous Lippmann-Schwinger equation to inverse scattering problems[END_REF].

However, the Factorization method approach, involved in [START_REF] Grisel | Determining the shape of defects in non-absorbing inhomogeneous media from far-field measurements[END_REF] and [START_REF] Bondarenko | The factorization method for inverse acoustic scattering in a layered medium[END_REF], provides a more straightforward formulation as an equivalence that is defined at each point. So, we propose here to use 1 this fast local information to reduce the computational effort in the complete refraction index reconstruction process.

More precisely, in this paper we are going to propose two hybridizations between Gauss-Newton and defect localization methods. We first consider the case where a known index has been locally modified. This could happen, for instance, from a local deterioration or a partially incorrect estimation of the actual index. In this case, a preliminary defect localization provides a geometrical information that we can use to choose which parameters have to be reconstructed. Then, the reconstruction can be performed straightforwardly on a reduced computational domain. The goal of this strategy is to avoid the useless reconstruction of parameters for which we have a suitable initial guess. In a second time, we investigate adaptive refinement. Here, defect localization is used to exhibit inaccurate regions in the current reconstruction. This local information allows us to refine the reconstruction mesh in these regions and resume the reconstruction to get a better precision while restraining the number of computed parameters. With this strategy, we thus aim to control the numerical effort by diminishing the number of discrete parameter values with respect to an uniform mesh of the whole parameter space.

This paper is structured as follows: In section 2, we specify the mathematical setting. We then introduce the direct problem in section 2.1, followed in section 2.2 by the description of the inverse medium problem we are interested in. The numerical method on which we will build our enhancements is then described in section 2.3. Afterwards, the defect localization and its applications are presented in section 3. We show how to reduce the reconstruction domain in section 3.1, and the adaptive refinement process is detailed in section 3.2. Finally, we numerically illustrate the sequence of both applications, and furthermore on a non-trivial example, in section 4. We end the the paper by concluding remarks in section 5.

Presentation of the problem

We start by introducing the direct problem and the inverse medium problem, followed by its numerical treatment.

The direct problem

We consider time-harmonic acoustic waves, with a fixed wave number k, modeled by the Helmholtz equation [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]. Inhomogeneous media are then represented by an acoustic refraction index, denoted by n ∈ L ∞ (R d ). So, the total field, denoted by

u n ∈ L 2 loc (R d ), is assumed to satisfy ∆u n + k 2 n(x)u n = 0, x ∈ R d , (1)
where d is the problem's dimension (d = 2 or 3). We consider compactly supported inhomogeneities, and we denote by D the support of n(x) -1. We also denote by u i ∈ L 2 loc (R d ) an incoming wave satisfying [START_REF] Scott | A review of non-destructive testing of composite materials[END_REF] with n(x) = 1 . The total field is then the sum of this incoming wave and the wave scattered by the inhomogeneous medium, denoted by u s ∈ L 2 loc (R d ):

u n := u s + u i , (2)
where the scattered wave is assumed to satisfy the Sommerfeld radiation condition

∂ r u s = iku s + O |x| -d-1 2 . (3)
Then, the linear system Besides, the outgoing part of a wave has an asymptotic behavior called the far field pattern, denoted by u ∞ n ∈ C ∞ (Γ m ), and given by the Atkinson expansion [START_REF] Kirsch | An introduction to the mathematical theory of inverse problems[END_REF]Theorem 6.11]

u n (x) := u i (x) + γ e ik|x| |x| d-1 2 u ∞ n ( x) + O |x| -d-1 2 , x := x |x| ∈ Γ m , (4)
where Γ m denotes the set of measurement directions as a subset of the unit sphere S d-1 (see figure 1), and where γ only depends on the dimension and is defined by

γ := e iπ/4 √ 8πk if d=2, 1 4π if d=3.
Furthermore, for practical reasons, we will mainly consider scattered waves having a plane wave source. These plane waves are defined by

u i ( θ, x) := exp(ik θ • x),
where θ is a unitary vector in the set of incidence directions, denoted by Γ e as shown in Figure 1. We then denote the total field with a plane wave source of incoming direction θ by

u n ( θ, x), θ ∈ Γ e , x ∈ R d .
Lastly, the corresponding far-field pattern in the measurement direction x ∈ Γ m will be denoted by

u ∞ n ( θ, x), θ ∈ Γ e , x ∈ Γ m .

The inverse medium problem

We are interested in the reconstruction, from far-field data, of an inhomogeneous refraction index that will be denoted by n ⋆ ∈ L ∞ (D) throughout this paper. Note that all considered indices will be implicitly extended by 1 outside D, the compact support of n ⋆ -1.

Most of the methods used to solve this reconstruction problem are valid for domains that are (subsets of) Hilbert spaces. Thus, we consider an index-to-far-field mapping F :

L ∞ (D) → L 2 (Γ e × Γ m ) defined by F(n) := u ∞ n .
Besides, data are generally perturbed by noise or measurement errors. So, we assume that we only have access to

u ε ∈ L 2 (Γ e × Γ m ), the perturbed version of u ∞ n ⋆ satisfying u ε -u ∞ n ⋆ L 2 (Γe×Γm) ε u ∞ n ⋆ L 2 (Γe×Γm) . (5)
The usual approach to this problem is then to find n by minimizing the difference between F(n) and u ε . Thus, we define the data misfit by

J(n) := c 1 F(n) -u ε 2 L 2 (Γe×Γm) ,
where c 1 denotes a normalization constant (e.g.

c 1 = u ε -2 L 2 (Γe×Γm)
). Even so, this problem is not continuous, as is shown by the following proposition. So, even a small perturbation ε can lead to a minimizer very far from n ⋆ . Proposition 2.1. The non-linear problem consisting in "finding n ε minimizing J" is ill-posed in the sense of Hadamard.

Proof. The mapping F is compact, and thus cannot have a continuous inverse.

Indeed, it has been shown in [14, Proposition 2.1.14] that the mapping n → u n is continuous from L ∞ (D) into H 2 (D); and as such, from L ∞ (D) into L 2 (D). As a consequence, the same property holds for the mapping f : n → (n -1)u n .

Moreover, the asymptotic behavior of the Lippmann-Schwinger equation yields the following relationship [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]Chap. 8.4]:

F(n)( θ, x) = z∈D e -ik x•z k 2 (n(z) -1)u n ( θ, z), θ ∈ Γ e , x ∈ Γ m . (6)
It is well known that the linear operator f → z∈D e -ik x•z k 2 f (z) is compact. Hence, the non-linear operator F is the combination of a compact linear operator with a continuous mapping. Therefore, it is compact itself.

Iterative approximation by a piecewise constant index

As stated in the introduction, we try to recover the unknown index n ⋆ with help of piecewise constant functions. Hence, the indices will numerically be represented by N complex parameters (η i ) i=1...N associated to the same number of zones (Z i ) i=1...N , so n(x) = i=1...N η i 1 Zi (x), where 1 Zi (x) is the characteristic function of Z i and i=1...N Z i = D. Each zone is thus a set of connected elements in the underlying mesh used to solve the Helmholtz equation. Moreover, to avoid any possibility of inverse crime, the reconstruction will be led on a second mesh that is different from the one used to generate the data u ε . As a consequence, the zones associated to the unknown parameters will intersect the discontinuities of n ⋆ . It is thus strictly impossible to reconstruct n ⋆ exactly. This is illustrated in Figure 2, showing two 2D meshes that will be used in our numerical simulations.

In this setting, we use the well-known Gauss-Newton method applied to the cost function J with a standard Tikhonov regularization term [START_REF] Mojabi | Overview and classification of some regularization techniques for the Gauss-Newton inversion method applied to inverse scattering problems[END_REF]:

J(n) := c 1 F(n) -u ε 2 L 2 (Γe×Γm) + c 2 n -n 0 2 L 2 (D) ,
where c 2 > 0 is a regularization parameter and n 0 ∈ L ∞ (D) is an initial guess. The choice of this regularization parameter parameter is discussed in a large number of papers, see for example [START_REF] Farquharson | A comparison of automatic techniques for estimating the regularization parameter in non-linear inverse problems[END_REF][START_REF] Bazán | An improved fixed-point algorithm for determining a Tikhonov regularization parameter[END_REF] and references therein. Empirically, it seems that a few percent of the fidelity term c 1 F(n)u ε 2 are a decent initial guess for c 2 . Besides, assumptions on n 0 and c 2 for the convergence of this method are discussed in [START_REF] Bakushinsky | The problem of the convergence of the iteratively regularized Gauss-Newton method[END_REF][START_REF] Blaschke | On convergence rates for the iteratively regularized Gauss-Newton method[END_REF]. Hence, the index n ⋆ we are looking for is approximated by a sequence (n p ) p∈N of indices, defined iteratively through Algorithm 1.

The gradient of the cost function, required for the Gauss-Newton method, has the following integral representation.

DF(n p ) ⋆ DF(n p ) + c 2 c 1 id (n p+1 -n 0 ) = -DF(n p ) ⋆ F(n p ) -u ε -DF(n p )(n p -n 0 ) ,
where id is the identity matrix, and DF(n p ) ⋆ is the hermitian adjoint of the matrix DF(n p );

Lemma 2.2. The mapping F is differentiable and the differential DF evaluated at n ∈ L ∞ (D) and applied to the direction dn ∈ L ∞ (D) has the following integral representation

DF(n) dn : ( θ, x) → z∈D k 2 u n (-x, z)u n ( θ, z) dn(z) dz, θ ∈ Γ e , x ∈ Γ m . (7)
Proof. Expansion (4) shows that u ∞ n ( θ, •) is linear with respect to the scattered field (u n -u i )( θ, •). Furthermore, It has been shown in [START_REF] Ségui | Sur un problème inverse en diffraction d'ondes : identification de permittivités complexes d'un matériau à partir de données du champ proche[END_REF]Proposition 4.3.1] that the scattered field is differentiable with respect to n and that the differential of the index-to-scattered-field mapping evaluated at

n ∈ L ∞ (D), applied to dn ∈ L ∞ (D), is the function v s ( θ, •) ∈ L 2 loc (R d ) satisfying ∆ x + k 2 n(x) v s ( θ, x) = -k 2 u n ( θ, x) dn(x), x ∈ R d , (8)
and the Sommerfeld radiation condition [START_REF] Mojabi | Microwave biomedical imaging using the multiplicative regularized gauss-newton inversion[END_REF]. Note that, contrarily to n, the direction dn is extended by 0 outside D. Thus, F is differentiable, and its differential is defined on

C ∞ (Γ e × Γ m ) by DF(n) dn = v ∞ .
Now, let us denote by Φ n (z, x) the Green function of the Helmholtz equation ( 1). Multiplying ( 8) by Φ n (z, x), integrating over D, and using Green's formula, yields

v s ( θ, x) = z∈D k 2 Φ n (z, x)u n ( θ, z) dn(z) dz, x ∈ R d .
The asymptotic behaviour is then given by

v ∞ ( θ, x) = z∈D k 2 Φ ∞ n (z, x)u n ( θ, z) dn(z) dz, x ∈ S d-1 .
Finally, representation ( 7) is obtained by applying the mixed reciprocity principle: [20, equation (3.66)]).

Φ ∞ n (z, x) = u n (-x, z) (see
Remark 2.3. As described in [START_REF] Bakushinsky | The problem of the convergence of the iteratively regularized Gauss-Newton method[END_REF][START_REF] Blaschke | On convergence rates for the iteratively regularized Gauss-Newton method[END_REF], to ensure the convergence of the regularized non-linear Gauss-Newton method (Algorithm 1), one also requires the differential of F to be Lipschitz continuous, which is actually the case. Indeed, the reference [14, Proposition 4.3.1], invoked in the proof of Lemma 2.2, gives the twice continuous differentiability of the scattered field with respect to n, thus implying that F is twice continuously differentiable too.

Numerical example

Set-up To illustrate our reconstruction schemes in R 2 , we consider a disc D of radius 1 centered at the origin. The embedded perturbation Ω is then chosen as another disc centered at (0.3, 0.3), and of radius 0.3. The mesh corresponding to this geometry can be seen in Figure 2a. The (perturbed) index we are looking for is set to n ⋆ := 1.3 in D \ Ω and n ⋆ := 1.6 in Ω whereas the initial guess, corresponding to the last known state, is

n 0 := 1.3 in D.
The Gaus-Newton method is performed with the regularization parameter c 2 := 10 -2 (and

c 1 = u ε -2
L 2 (Γe×Γm) , as previously denoted). Also, solutions to the Helmholtz equation are computed via a P 1 finite element method and Cartesian Perfectly Matched Layers (PML) [START_REF] Dah | Sur la modélisation de milieux fictifs absorbants de type couches de Bérenger[END_REF]. Lastly, the corresponding far-fields are evaluated through the representation formula [START_REF] Grisel | Determining the shape of defects in non-absorbing inhomogeneous media from far-field measurements[END_REF]. For all these examples, the wave number is set to k = 5, and the angles corresponding to incoming/measurement directions are equally distributed over [0, 2π]. Results An example can be seen in Figure 3 with a reconstruction mesh of 2672 triangles (see Figure 2b) divided into N = 10, N = 27, N = 75 and N = 2672 zones. More precisely, the final relative error

e p End := n p End -n ⋆ L 2 (D) / n ⋆ L 2 (D)
is synthesized as a function of the number of zones N in Figure 4. Moreover, for comparison purpose, we list in Table 1 the exact final relative error obtained in several configurations. Besides, we see in this table that the relative error is of order 10 -2 , so we choose the stopping criterion n p -n p-1 2 /(1+ n p-1 2 ) < ǫ = 10 -4 for all our reconstructions. In all these test cases, this bound was reached after four iterations.

Comments The low error that can be seen in Figure 4 for N = 19 is a particular case related to the considered test case. Indeed, it just happens that this choice of zones provides a natural match to our simple geometry, yielding a reconstruction that is better than expected. In the presented piecewise constant iterative approximation, the possible precision is directly linked to the number of basis functions N which, in turn, is linked to the computational effort.

p En d ) in % ε = 5% ε = 2% ε = 1%
In the lack of some geometrical informations, all parameters are equally treated and updated at each iteration. However, this can generate more effort than is really needed, and we address two cases where these unnecessary efforts can be reduced. For the second case, we are concerned in focusing on the most inexact constants during the reconstruction. Indeed, to obtain a precise identification, the reconstruction mesh has to be refined in the zones intersected by the discontinuities of n ⋆ . However, if n ⋆ is constant in large areas, refining the reconstruction mesh everywhere only raises the computational effort for a relatively small precision increment.

To address these aspects of the reconstruction, the useful information in both cases would thus be the localization of the nearly exact constants. Or, to put it otherwise, we look for the localization of the difference between the reconstruction at hand and the exact (unknown) index. We call this defect localization with respect to an inhomogeneous background reference. Of course, to enhance the complete reconstruction, access to this specific information should be fast. To this end, it has been shown in [START_REF] Grisel | Determining the shape of defects in non-absorbing inhomogeneous media from far-field measurements[END_REF] that the Factorization method for a constant background, as presented in [START_REF] Kirsch | The factorization method for inverse problems[END_REF], can be extended to an inhomogeneous background with help of a modified measurement operator. . Lastly, assume that n and n ⋆ are real valued, and that either (n -n ⋆ ) or (n ⋆ -n) is locally bounded from below in Ω := support(n -n ⋆ ).

Then, for each z ∈ R d , we have n(z) = n ⋆ (z) if, and only if,

S {n,n ⋆ } (z) :=    j u n (•, z), ψ j L 2 (S d-1 ) 2 σ j    -1 > 0,
where (σ j , ψ j ) is an eigensystem of W # .

Remark 3.2. This localization result calls for a number of remarks.

1. Theorem 3.1 requires full bi-static data (Γ m = Γ e = S d-1 ) and real-valued indices. However, we also recall the conjecture, stated in [6, Remark 6.2]: To build the localization function S, the eigensystem of W # , denoted by (σ j , ψ j ), could be replaced by a right-singular system of (F n ⋆ -F n ). The main benefit is the possibility of considering Γ m = Γ e = S d-1 and complex valued indices. Moreover, we a priori need to truncate the sum to avoid overvalued terms, due to perturbations on the smallest eigenvalues. However, in our experience, we have noticed that the computations are stable. So, no truncation is performed for the numerical applications of this paper and all the eigenvalues of the measurement operator's matrix representation are considered.

2. Furthermore, numerical examples in [START_REF] Grisel | Determining the shape of defects in non-absorbing inhomogeneous media from far-field measurements[END_REF] show that this localization is effective for defects bigger than (approximately) one over six of the wavelength. Besides, in order to get satisfactory results in the successive resolutions of the Helmholtz equation, we have set the reconstruction mesh size to be about one over twenty of the wavelength. Thus, we will only consider defects that cover at least four connected mesh elements.

3. Moreover, the examples shown in [START_REF] Grisel | Determining the shape of defects in non-absorbing inhomogeneous media from far-field measurements[END_REF] exhibit that defects can be localized even when the surrounding background is not precisely known. Practically, low amplitude inaccuracies with respect to the exact index do not seem to interfere with the localization of the contrasting defects. Thus, geometrical information gained through the defect localization presented here is expected to focus on the most "defective" zones.

4. Finally, it is to be noted that the added computational cost of this localization function within the Gauss-Newton process will be negligible. Indeed, the functions u n (•, z) will already have been computed to evaluate the differential of F, as can be seen on the integral representation (7).

Selective reconstruction

We here consider the case where the initial guess n 0 is exact, except for some perturbation whose support will be denoted by Ω. Thus, we propose to perform a preliminary selection of the parameters, to reconstruct only the perturbed ones. The selection is performed by considering only the parameters associated to zones where the maximal value of the (normalized) defect localization function S {n0,n ⋆ } / max D S {n0,n ⋆ } is above some threshold T . The whole index n ⋆ is then reconstructed by updating those parameters only. This leads to a reconstruction, described in Algorithm 2, using a number of parameters N Sel that should be significantly less than N .

Algorithm 2: Selective reconstruction Input:

n 0 ∈ L 2 (D) 1 S i ← max Zi S {n0,n ⋆ } (x);
2 Ω T ← the set of zones for which S i > T max S i ; 3 n p End ← Algorithm 1(n 0 Ω T ) (all indices are extended by n 0 outside Ω T );

Output: n p End

Numerical example

Set-up In the framework of section 2.3, we here consider the smallest possible zones, that is one parameter for each triangle of the reconstruction mesh. Figure 5 shows which zones are selected with three threshold values T = 10%, T = 20% and T = 30%.

Results

We can see that a threshold of T = 10% yields an accurate selection of the perturbation, and thus provides a satisfactory reconstruction with only N Sel = 323 selected parameters. Thus, we end up with significantly less parameters than the 2672 we have initially considered.

Comments We can also see in Table 2 that the relative error can be lower than what was obtained through a full Gauss-Newton reconstruction over a set of various configurations. This is a consequence of the fact that all the parameters outside the perturbation are equal to the exact value, while they can be miscalculated in the full reconstruction. Identifying the unperturbed parameters can thus clearly enhance the reconstruction. As previously, the stopping criterion was reached after four iterations in all cases. However, a threshold of T = 20% seems too high, as the 181 selected zones do not completely cover the perturbation's support, resulting in a slightly flawed reconstruction. More precisely, the relative error obtained as a function of T with 30 × 30 data can be seen in Figure 6. Clearly, there is an optimal value for T around 10% when the noise ratio is kept low.

Besides, with more noise (5%), we see in Figure 6 that the optimal T is shifted towards 20%. Furthermore, we see that a good estimation of this threshold becomes even more important when the noise level grows. This brings up the problem of how to select a correct threshold, taking at least the measurements noise and the amount of data into account. Unfortunately, for the moment, we do not have a realistic indicator to tell if the selected threshold is acceptable.

Adaptive refinement

As stated in section 2.3, we use a reconstruction mesh that is different from the one used to generate the data. Hence, the supports of the basis functions used in the reconstruction will not follow the geometry of n ⋆ , especially with a low number N of basis functions. Thus, we propose to iteratively refine the reconstruction mesh with help of the previously introduced defect localization, in order to provide a satisfying approximation of the unknown index with a small number of parameters. The refinement outline is presented in Algorithm 3.

The number of 16 mesh elements is taken so that, after the splitting, each zone has still more than four mesh elements, which is the lower limit for defects to be relevant, as specified in Remark 3.2.

Numerical example

Set-up We illustrate our adaptive refinement in Figure 7 in the same conditions as in section 3.1. Selection Threshold (T ) in % Final relative error (e

p En d ) in % ε = 5% ε = 2% ε = 1%
n 0 ∈ L ∞ (D) 1 p ← 0; 2 repeat 3 S i ← max Zi S {np,n ⋆ } (x); 4 I ← {i such that Z i contains

Results

The steps 3 (defect localization) and 9 (reconstruction on the refined set) of Algorithm 3 are illustrated alternately in Figures 7a-7f, and it can be seen how the reconstruction focuses on the support of the contrasting perturbation. Figure 7g represents the values of n 59 , which is obtained with N = 76 basis functions chosen during 25 successive adaptive refinements. Also, the relative error e p , obtained in step 10 of the algorithm, is plotted in Figure 7h as a function of p.

Comments First, it can be noted that each refinement adds 3 parameters to be reconstructed and that each call to Algorithm 1 generates about four iterations (see Tables 12). So, the number of iterations is comparable to the number of parameters. Then, comparing with the results obtained when using basis functions that are placed randomly, summarized in Table 1 or in Figure 3f, we can see lower reconstruction errors when using our guided adaptive refinement. In this example, our results are even comparable to the complete reconstruction (Algorithm 1) performed with 20 times more parameters. We thus meet our goal, which is a satisfactory reconstruction with a limited number of well-chosen basis functions.

Moreover, as we can see in Table 3, the sensitivities to noise or data amount in this example are similar to what we observe in section 3.1.

Finally, it cannot be overlooked that the number of total iterations p end is now quite high, since each loop in Algorithm 3 computes an iterative reconstruction. However each of those reconstructions is conducted on a very small number of parameters. So, the integral representation (7) for the differential is not cost effective in this case and a suitably tuned evaluation of n p+1 in Algorithm 1 might thus be able to balance the higher number of iterations.

Combining both strategies

The selective reconstruction is presented in section 3.1 as a preliminary step to the reconstruction. Furthermore, the adaptive refinement described in section 3.2 enhances the actual reconstruction step. So, adaptive refinement and selective reconstruction can be used one after the other. This extension of Algorithm 2 is described in Algorithm 4.

Algorithm 4: Selective reconstruction followed by adaptive refinement

Input: n 0 ∈ L ∞ (D) 1 S i ← max Zi S {n0,n ⋆ } (x); 2 Ω T ← the set of zones on which S i > T max(S i ); 3 n p End ← Algorithm 3(n 0 Ω T ) (all indices are extended by n 0 outside Ω T ); Output: n p End
Note that the number of parameters selected in step 1 of this algorithm is not directly used in the adaptively refined reconstruction (step 3). Indeed, the iterative refinement described in Algorithm 3 starts the reconstruction with only one zone. More precisely, the information retained from the selection step is the shape of the perturbation. Note that the accuracy of this selection is important: This is what allows the adaptive refinement to focus on the reconstruction of the perturbation's inner geometry, instead of focusing on the contrast between the perturbation and the background.

Numerical example 1

Set-up As in section 3, we illustrate Algorithm 4 with the selection thresholds T = 10%, T = 20% and T = 30%. The respective selected mesh elements can be seen in Figures 5a-5e.

Results Figures 8a-8c show the reconstructions after 2, 4 and final adaptive refinement loops with a threshold T = 10%. As expected through the previous results, the reconstruction is very good. In fact, the exact values listed in Table 4 show that this reconstruction reaches an accuracy comparable to the one obtained through the initial selective reconstruction; the latter requiring 10 times more basis functions. As in section 3.2, and for the same reasons, the number of parameters for each adaptively refined reconstruction is comparable to the number of iterations.

Similarly to the examples presented in section 3.1, T 20% also provides a too small selection, leading to a flawed reconstruction. For example, it can be seen in Figures 8d-8i that the reconstruction tends to a crown shape. So, the constraint induced by this too small selection seems to create false local minima, altering he whole convergence process, as we can see in Figure 8j. The corresponding relative error values are detailed in Table 4.

Comments Since the selection is performed before the adaptive refinement, the choice of the threshold T still has a large influence in the final result, even in the case of over-selection. However, results in terms of accuracy remain close to the reference listed in Table 1 while involving only 0.6% to 2% of the total number of elements used in the full Gauss-Newton reconstruction. Set-up As a last example, we now consider a more elaborate and complex valued unknown index n ⋆ , shown in Figure 9. Besides, we also make this reconstruction more challenging by reducing the measurements aperture. Incoming directions are still taken in [0, 2π], but there will be five less, and measurement directions are now taken in [0, 3 2 π]. In this situation, the 1 cannot be defined. So, we consider the technical modification, recalled in Remark 3.2, that is conjectured to cover this case. Furthermore, we assume that n ⋆ was known before the central perturbation. So, we consider the initial guess n 0 shown in Figure 10. Finally, to remain in the previously defined context, we present the results of Algorithm 4 applied to this new geometry with the same selection thresholds T = 10%, T = 20% and T = 30%.

Results

The reference reconstructions obtained with the usual Gauss-Newton reconstruction (Algorithm 1) in the special case of 30 × 25 data and 2% noise are synthesized in Figure 11.

We then present in Figure 12 the selected zones and the resulting reconstruction corresponding to each selection threshold. In this case, T = 20% now seems to be the best threshold value, and T = 30% is still too high. Indeed, in the lines of the previous numerical example depicted in Figure 8, T = 30% seems again to induce a false local minima, resulting in an altered convergence sequence. This is confirmed in Figure 12j, where we can also see that, even though T = 10% allows to reach a satisfying precision, it requires much more refinements to do so than with T = 20%.

The results obtained in section 3 are thus reinforced by this example, exhibiting reconstructions comparable in precision to the full Gauss-Newton reconstruction, but with a much lower number of parameters. Comments Note that with this less trivial test case, the borders of the supports of the basis functions for the reconstruction do not correspond to the discontinuities of the exact index. This ensures that we avoid this particular case, which was mentioned in the concluding comments of Section 2.3.

Conclusion

We have used a defect localization method to propose two ways of reducing the number of parameters in the reconstruction of an unknown refraction index. The first method is set in the context of defects identification and uses their localization to reconstruct only the useful parameters of the whole index. The second method is an adaptive refinement, based on defect localization to iteratively reconstruct a better approximation with a limited number of parameters. We have obtained good numerical results with both methods. The reconstruction could however be further enhanced by two automations: some automatic choice of the threshold for the defect localization function and some automatic selection of the regularization parameter. The second issue has been reviewed for example in [START_REF] Farquharson | A comparison of automatic techniques for estimating the regularization parameter in non-linear inverse problems[END_REF][START_REF] Bazán | An improved fixed-point algorithm for determining a Tikhonov regularization parameter[END_REF] and is claimed to be less critical when using a so-called Multiplicative Regularization described in [START_REF] Van Den Berg | Extended contrast source inversion[END_REF]. However, for now, we have not been able to further enhance our results with these techniques. Finally, convergence of the coupled process, hybridizing the Gauss-Newton and Factorization methods, presented in this paper remains to be investigated.

Figure 1 :

 1 Figure 1: General setting and notations.

Figure 2 : 3

 23 Figure 2: Test case geometry

  Evolution of the relative error

Figure 3 :

 3 Figure 3: Gauss-Newton reconstruction with 30 × 30 data and ε = 2% noise

  Numb er of parameters (N ) Final relative error (e

Figure 4 :

 4 Figure 4: Final relative error with 30 × 30 data and different noise levels ε

Theorem 3 . 1 .

 31 [START_REF] Grisel | Determining the shape of defects in non-absorbing inhomogeneous media from far-field measurements[END_REF] Theorem 6.1] Assume that Γ m = Γ e = S d-1 . Then, define a measurement operatorW := Id + 2ik |γ| 2 F n (F n ⋆ -F n ) , where F n : L 2 (S d-1 ) → L 2 (S d-1 )denotes the classical far-field operator, defined by F n g( x) := g, u ∞ n (•, x) L 2 (S d-1 ) . Next, we define the positive self-adjoint operator W # by W # := |W + W ⋆ | + |W -W ⋆ |, where the notation |•| applied to an operator L stands for |L| := (L ⋆ L) 1 2

Figure 5 :

 5 Figure 5: Selective reconstruction with 30 × 30 data and ε = 2% noise

Figure 6 :

 6 Figure 6: Influence of the threshold T with 30 × 30 data and different noise levels ε

  Relative error (ep) in %(h) Evolution of the relative error

Figure 7 :

 7 Figure 7: Adaptive refinement (Algorithm 3) with 30 × 30 data and ε = 2% noise

  Evolution of the relative error

Figure 8 :

 8 Figure 8: Selective reconstruction followed by adaptive refinement, 30 × 30 data and ε = 2% noise 17

Figure 9 : 6 (b) Im n 0 Figure 10 :

 96010 Figure 9: Exact index n ⋆

Figure 11 :

 11 Figure 11: Relative error for the usual Gauss-Newton method, with 30 × 25 data and different noise levels ε

Figure 12 :

 12 Figure 12: Selective reconstruction chained with adaptive refinement for the more elaborate example with 30 × 30 data and ε = 2% noise

Table 1 :

 1 1. For the first case, we consider a bounded perturbation in a known initial state n 0 . So, most of the values of the index have not changed and should not be reconstructed. Gauss Newton reconstruction 2.

			15 × 15 data 30 × 30 data 60 × 60 data
	N	ε	e p End	e p End	e p End
		5%	5.9%	5.9%	5.9%
	10	2%	5.9%	5.9%	5.8%
		1%	5.9%	5.8%	5.8%
		5%	4.9%	4.9%	4.9%
	27	2%	4.8%	4.8%	4.8%
		1%	4.9%	4.8%	4.8%
		5%	5.4%	5.0%	4.4%
	75	2%	3.9%	3.7%	3.7%
		1%	3.6%	3.6%	3.6%
		5%	5.3%	4.5%	3.9%
	2672	2%	3.5%	3.3%	3.1%
		1%	3.1%	3.0%	2.9%

Table 2 :

 2 Sel e p End N Sel e p End N Sel e p End Selective reconstruction

		15 × 15 data 30 × 30 data 60 × 60 data
	T N 10% ε 5% 874 2% 354	4.0% 2.4%	739 323	3.3% 2.3%	633 360	2.8% 2.3%
		1% 305	2.3%	282	2.4%	296	2.3%
		5% 321	2.7%	282	2.6%	268	2.8%
	20%	2% 196	3.5%	181	3.7%	203	3.3%
		1% 172	4.1%	162	4.3%	171	4.0%
		5% 204	3.3%	181	3.7%	178	3.9%
	30%	2% 134	5.4%	125	5.7%	136	5.3%
		1% 120	5.8%	112	5.9%	115	5.8%

  more than 16 mesh elements}; Split ← i such that S i Split = max i∈I S i ; Update the set of zones by splitting Z i Split into four sub-zones; Update the set of parameters accordingly by duplicating η i Split three times; End ; 11 until N > N max or each Z i contains less than 16 mesh elements; Output: n p End

5 i 6 7 8 N ← N + 3; 9 n p+p End ← Algorithm 1(n p ); 10 p ← p + p

Table 3 :

 3 Adaptive refinement

		15 × 15 data 30 × 30 data 60 × 60 data
	ε	N	e p End	N	e p End	N	e p End
	5% 76	4.9%	76	6.0%	76	5.0%
	2% 76	4.0%	76	4.8%	76	3.3%
	1% 76	3.8%	76	4.4%	76	3.6%

Table 4 :

 4 Selective reconstruction chained with iterative refinement

			15 × 15 data 30 × 30 data 60 × 60 data
	T	ε	N	e p End	N	e p End	N	e p End
		5% 60	3.9%	57	3.1%	61	3.1%
	10%	2% 52	2.6%	52	2.4%	55	2.5%
		1% 49	2.6%	52	2.5%	52	2.5%
		5% 52	3.3%	52	2.8%	46	3.0%
	20%	2% 22	3.8%	19	3.7%	22	3.3%
		1% 19	4.0%	19	3.9%	16	4.1%
		5% 16	4.1%	19	3.6%	16	4.2%
	30%	2% 16	5.3%	16	5.6%	16	5.3%
		1% 15	5.5%	16	5.6%	16	5.6%
	Numerical example 2

p ← p + 1;

until n p -n p-1 2 /(1 + n p-1 2 ) < ǫ; Output: n p End

Acknowledgments

Support for some of the authors of this work was provided by the FRAE (Fondation de Recherche pour l'Aéronautique et l'Espace, http://www.fnrae.org/), research project IPPON.