A general method to reconstruct strong gravitational lenses based on the singular perturbative approach

The number of gravitational arcs systems detected is increasing quickly and should even increase at a faster rate in the near future. This wealth of new gravitational arcs requires the development of a purely automated method to reconstruct the lens and source. A general reconstruction method based on the singular perturbative approach is proposed in this paper. This method generates a lens and source reconstruction directly from the gravitational arc image. The method is fully automated and works in two steps. The first step is to generate a guess solution based on the circular solution in the singular perturbative approach. The second step is to break the sign degeneracy and to refine the solution by using a general source model. The refinement of the solution is conducted step by step to avoid the source-lens degeneracy issue. One important asset of this automated method is that the lens solution is written in universal terms which allows the computation of statistics. Considering the large number of lenses which should be available in the near future this ability to compute un-biased statistics is an important asset.


INTRODUCTION
There are basically two main issues with the automated reconstruction of gravitational arcs: (1) what model to choose for a given gravitational arc ?(2) Given the large parameter space for the models, how to identify some first guess for the solution ?We will show in this paper that the singular perturbative approach offers an efficient solution to these two issues.The paper starts with an introduction to the perturbative singular approach, and continues with a discussion on the proper referential for the reconstruction of the lens.In this referential a first guess for the solution is constructed based on the circular source model in the perturbative approach (see Alard (2007) Eq. ( 12)).Once this first guess is reconstructed the next step will be to use a general chi-square minimization method to refine the former guess and reach the optimal χ 2 minimum.The degeneracy of the first guess is naturally broken by constructing a refined non circularly symmetric source model.The reconstruction is completed by increasing the potential Fourier order to the required minimal order.The first part of the paper will recall the basics of the singular perturbative approach.This will be followed by a description of the methods and associated equations.The presentation will be illustrated with a number of practical applications by reconstructing simulated gravitational arcs images.

BASICS OF THE SINGULAR PERTURBATIVE APPROACH IN GRAVITATIONAL LENSING.
We define the local coordinates r and source coordinates r S , in a lens potential φ the lens equation reads, The basic idea of the singular perturbative approach is to consider that the un-perturbed situation corresponds to the image of a point at the center of a circularly symmetric potential.
Due to its symmetry the image of a point in a circularly symmetric potential is a circle corresponding to the Einstein circle (R E ).For convenience we re-normalize the coordinate system to obtain R E = 1.As consequence the unperturbed situation is an infinity of points located on the Einstein circle.The choice of this un-perturbed situation allow us to deal with the high non-linearity of strong gravitational lenses in the angular direction.In effect at any angular position θ an un-perturbed point is present.On the other hand, the perturbation in the radial direction is generally small for gravitational arcs, which allow us to work in the vicinity of the Einstein circle.It is also assumed that strong gravitational arcs are due to small perturbation of a circular potential.Consequently in this framework we have, φ (r, θ ) = φ 0 (r) + εψ(r, θ ) r = 1 + εdr r S = εr s (2) MNRAS 000, 1-?? (2023) In Eq. ( 2) ε is a small number.By introducing Eq. ( 2) in the lens equation Eq. ( 1) the perturbative lens equation is obtained (see Alard (2007) Eq. ( 8) and ( 7)), (3) It is useful to include the source impact parameter r 0 in the definition of the perturbative fields f 1 and d f 0 dθ by defining the new fields f1 and d f0 dθ , In the continuation of this paper we will assume the impact parameters are included in the field and will not use the fi symbol to lighten the notation and will use directly the f i notation.

2.1
The solution for a circular source in the perturbative approach.
An interesting application is to consider the image of a round source contour with radius R 0 in the singular perturbative theory.This model is simple with a direct analytical solution for the images.In effect in this case solving Eq. ( 3) is straightforward and reduces to a second order equation (see Alard (2007), Eq. 12).
It is interesting to note that in this circular source model the two perturbative fields have simple physical meanings.The field f 1 represents the mean position of the image for each angular position, while the field d f 0 dθ is directly related to the image width.

Relation to the multipole expansion
A final and interesting point about the singular perturbative expansion is that the perturbative fields are related to the multipole expansion on the Einstein circle R E (here re-scale at R E = 1).The multipole expansion of the potential at R E = 1 reads: The coefficients, a n (r), b n (r), c n (r), d n (r) of the potential multipole expansion (Eq.6) are related to the surface density Σ by the following equations, This multipole expansion of Eq's ( 6) and ( 7) is related to the singular perturbative expansion by Eq. ( 8) (see Alard (2009) for the details of the calculations.) It is interesting to note that Eq. ( 8) implies a direct relation between the Fourier expansion of the perturbative fields and the multipole expansion of the potential on the Einstein circle.
This also implies that in this theory the inner and outer contribution to the lens potential can be separated.An example of a reconstruction with a demonstrative separation of the inner and outer potential contribution can be found in Alard (2010)).

THE METHOD
3.1 Estimating the center of the perturbative reconstruction.
A basic problem in the perturbative reconstruction of the lens is that the actual center of the unperturbed circular distribution is unknown.In most case the lens itself can be observed and one could use the center of the light distribution in the lens as the center for the reconstruction.However even if this assumption seems reasonable it might be wrong in some cases.It is more interesting to develop a general method without making assumptions and find the center as a parameter of the reconstruction.We will start by considering that we are working in a referential which is slightly shifted form the optimal (true) referential.
The effect of the small shift between the two referential on the perturbative fields will be estimated using a perturbative expansion.
3.1.1Effect of a small shift on the perturbative fields.
Let's define the offset from the optimally centered referential, δ 0 .We will assume that, The coordinates in the optimal referential and offset referential are respectively r and u.The two coordinates are related by, The lens equation (Eq.( 1)) is transformed using Eq. ( 10) leading to, By expanding Eq. ( 11) to first order in ε with the help of Eq. ( 23) we obtain, Considering that φ = φ 0 + εψ (see Eq. 2) and again expanding to first order in ε we obtain, According to Eq. (3) we have, The notation ψ represent the field ψ once the effect of the offset is taken into account.As a consequence for the offset fields f 1 and d f 0 dθ using Eq. ( 13) we obtain, Note that in Eq (14) we used ∇φ 1 = δ to obtain the radial and angular derivatives of φ 1 (see Eq. ( 13)).We also used the fact that dφ 0 dr r=1 which is the consequence of the critical condition for the unperturbed potential at the scale radius R E = 1.Here we used the variable κ 2 which was already used in the initial presentation of the singular perturbative theory, see Alard (2007).It is interesting to note that to the first order Eq. ( 14) indicates that the field dθ is unaffected by the shift δ 0 .The field f 1 itself is affected by the shift.

3.2
The choice of the referential.
In general the Fourier expansion of the perturbative fields at first order contain four independent terms.Let's define these first order terms as α and β for f 1 and d f 0 dθ .We have also to include the source impact parameter r 0 (see Eq. 4).The first order Fourier of a perturbative fields F is represented by the suffix We will consider that in our optimal referential these first order (α, β ) Fourier terms are not zero but are the smallest possible.If we consider a slight change in referential by introducing a small shift λ by using Eq.14 we obtain, It is interesting to note that Eq. ( 16) implies a fundamental degeneracy in the choice of the referential.In effect one could choose λ such that, Since in the optimal referential model |β | is small then by direct consequence of Eq. ( 17) the coefficients at first order for f 1 are small two.Thus by this small shift we found another referential where the first order terms are small.The basic idea behind this choice is that most dark matter halos are not too distant from an isothermal profile.Furthermore it is simple to show that for slightly non-isothermal halos the fields are conjugates for first order Fourier terms.

Expression of the perturbative fields for slightly non-isothermal potentials.
For a centered circular isothermal potential both perturbative fields are null.When the same potential is slightly off-centered and shifted by λ at the first order the potential will be, The corresponding perturbative fields can be estimated by using f 1 = dφ dr r=1 − 1 + r 0 u r and + r 0 u θ .Here r 0 is the source impact parameter (see Eq. 4).
If we consider a set of slightly of centered isothermal with corresponding weights potential Eq. ( 19) indicates that the contribution to the perturbative fields will be the weighted sum of of the of the impact parameters of each component.Thus for such slightly perturbed isothermal potentials we see that the perturbative fields will have conjugate expressions and depends only on two terms.As a consequence we will choose the conjugate field model in the optimal referential.This is a reasonable choice for dark matter halos in this degenerate problem.In any case as long as the shift from the optimal referential is small and of the order of the perturbation size, the specific choice of the referential should not matter and should not affect the final result.
3.3 Estimating a first guess for the solution.
The process of building the initial guess is de-composed in two steps, first an estimation of the Einstein ring size and center, and second an estimation of the first guess of perturbative fields based on a round source model.

Estimating the position and size of the Einstein ring.
The first step is make some estimation of the radius and center of the Einstein circle associated to the lens.The most direct method is to fit a circle to the arc system.To fit a circle to the data we will consider a model where the arc system is a slightly off-centered in the current coordinate system.As a first guess for the center of the circle we may consider the center of gravity the arc system or more directly the center of the image.The actual choice of the guess is not very important as the search for the center and radius of the circle converge quickly even for quite a distant initial guess.For a circle of radius R c off-centered by r c to the first order in r c R c the circle equation reads, The first step of our processing is to reduce the image of the arc to a set of positions defined by the radius and corresponding angle.The processing starts by choosing a threshold.In the calculations we will consider only the pixels with values above this threshold.In practice the angular domain is divided into a number of bins N θ and for each bin the mean radius R θ is estimated.Then the First order Fourier series defined by Eq. ( 20) is fitted to the set of positions (R θ , θ ).This fitting procedure estimates the parameters r 0 , and r c .Note that first we estimate the parameter α in Eq. ( 20) and then we infer r c .Once the estimation of the parameters has been performed the referential is shifted according to the estimation of r c .The the process is iterated, a new estimation is performed and the center is refined.
This fitting process is very efficient and converges quickly.The result is a first estimation of the Einstein radius and a first estimation of the lens center.The circle fitting procedure converges when the circle is centered in the referential system.When the method converges and the circle is centered in the referential system the first order Fourier terms in Eq. ( 20) are zero.For a circular source model the field f 1 represents the mean position of the images (see Eq. 5 and Sec.2.1).As a consequence for a simple round source model the first order Fourier terms of the fields f 1 have to be zero in the final referential after the circle fitting procedure has converged.In a system shifted by a vector δ 0 from the optimal referential Eq. ( 14) and Eq. ( 19) implies that, In Eq. ( 21) the degenerate term λ + r 0 was re-defined by introducing the r0 = λ + r 0 .Note that in the referential defined by the circle fitting procedure the first order terms have to be zero, which by using Eq. ( 21) is equivalent to, Our choice for the centering of the arcs is based on an isothermal halo conjugated model (See Sec.3.2.1).For an isothermal halo we have, κ 2 = 1, r0 = r 0 , Eq. ( 22) implies that δ 0 = r 0 which means that the shift δ 0 corresponds to the source impact parameter.As a consequence for an isothermal potential the circle centering procedure implies that the referential is centered on the source.As illustrated by the rotation curves the dark matter halos are not far from isothermal thus we the referential will be close to a source centered referential.An interesting point is that at first order the field d f0 dθ is un-affected by the shift and consequently allows a direct estimation of the pseudo impact parameter r0 (see Eq. ( 22)).

Estimating a first guess for the perturbative fields.
A first rough reconstruction of the perturbative fields in this lens centered initial referential will be performed.In the continuation the estimation of the angular bin mean positions and width will be all renormalized by the estimated Einstein radius.In the circular source model the field f 1 is simply the image mean position while the field d f 0 dθ is related to the image The reconstruction of the field f 1 is obtained directly by fitting a Fourier series to the mean position R θ of the binned data at each angular bin position where an image is available.
For the field d f 0 dθ the data to fit will be the image width ∆ (Eq.23) for each angular bin.An estimation of ∆ will be obtained by counting the number of pixel in each angular bin above the threshold.The number of pixels above the threshold is a direct estimation of the bin surface from which the effective width ∆ can be easily estimated.In practice the number of pixel above a threshold is an integer quantity and as a consequence may not be accurate for a small number of pixels.To improve the accuracy of this estimate the number of pixels just below the threshold and just above the threshold are estimated (see Eq. ( 26)).
An interpolation of these two nearby pixel counts provide an accurate estimation of the pixels number at the threshold for each angular bin C θ .
It is in principle possible to evaluate the field d f 0 dθ by relating the numerical data of Eq. ( 24) with equation Eq. ( 23).However in practice a direct evaluation of d f 0 dθ using Eq's ( 23) and ( 24) is problematic for two reasons.The first reason is that due to the square in Eq. ( 23) the sign of d f 0 dθ is unknown.Each time the function d f 0 dθ cross the zero line the sign may or may not change.The field d f 0 dθ is zero near the center of each image, thus if there are many images the uncertainty on the sign means that the number of potential solution to explore is large.And furthermore the second issue is that the management of the noise is very delicate due to the square root in Eq. ( 23).An optimal estimation of d f 0 dθ requires a fitting procedure based on a least-square model.The Fourier series are the simplest and most direct non parametric reconstruction of the perturbative fields since they are directly related to the multipole expansion on the Einstein circle (see Eq. ( 6)).To construct the initial guess it is important to reduce the complexity of the solution to a minimum and thus to consider an expansion in Fourier series of the fields of lower order.A typical number for the order of the Fourier N F expansion is N F = 2 or N F = 3 depending on the complexity of the solution.It would be very unusual and quite impractical to reach N F = 4 for the initial guess.Another point is the existence of a possible mass-sheet degeneracy related to the scaling coefficient κ 2 (see Eq. ( 3)).In order to cancel the effect of the re-scaling coefficient κ 2 it is useful to introduce the correlation coefficient which is scale free rather than a least-square estimate.
As a consequence the estimation of d f 0 dθ will be performed by maximizing the correlation coefficient between the two expression of ∆, the analytical expression in Eq. ( 23) and the numerical data of Eq. ( 24).The order of the Fourier series is low and the cost of computing the correlation coefficient for the array of binned data is low since in practice the number of bins is less than one hundred.As a consequence it is definitely practical to explore the parameter space to estimate a model for d f 0 dθ .The field d f 0 dθ has no zero order term since it is a derivative, thus at order 2 the number of coefficients is 4 and 6 at order 3.The parameter space is explored by generating a series of random numbers in a typical range.The range to investigate is by nature of a fraction of unity.In practice the range at order 1 is about 0.5, and decrease by a scale factor corresponding to the order at higher order.The assumption of decreasing order is justified by the fact that the power of the Fourier series at some higher order will be close to zero.Thus between the first order and this higher order a decrease is observed and we make the choice of the simplest model, a linear decrease.In practice this choice is notation very important since we start with a very broad range at order 1 which is much larger than practical cases.This model of linear decrease can thus be considered as an upper bound.At order 2 the exploration of the parameter space is performed using 10 6 set of random number.The corresponding computing time is only of about a few seconds on a basic PC for 75 angular bins.Considering the low computing cost the number of random sets can be easily increased if higher order Fourier solutions have to be explored.

Illustration using a numerical simulation.
The application of this method to a numerical simulation of a cusp configuration for an elliptical NFW potential is presented in

Refining the guess and converging to an optimal solution.
Once a a first approximate solution based on the circular source model has been obtained for the fields the problem is to refine this solution and find the optimal solution.In practice the approximate solution will be used as a starting guess for a full minimization in the multi-parameter space.However the first thing to consider in the search for an optimal solution is to define an optimal referential.An optimal referential should minimize the noncircular parts of the lens potential and as a consequence should be closely centered on the lens.However the circle fitting procedure implies that the first guess is defined in a source centered referential (see Sec. (3.3.1) and Eq's ( 21) and ( 22)).Consequently the first step should be to move from the source centered referential where we defined the initial guess to a lens centered referential.This task is facilitated by the fact that the field d f 0 dθ is unchanged by a small shift (see Eq. ( 14)).When expressed in a referential shifted by δ 0 the first order term of d f 0 dθ are precisely the pseudo-impact parameter r0 (see Eq. ( 21)).Thus from the guess solution for d f 0 dθ we can estimate the pseudo-impact parameter ¯r0 .A problem is that the guess solution is estimated from the circular solution of Eq. ( 5).In this equation the square of d f 0 dθ appear which mean that a degeneracy on sign of d f 0 dθ is present.This means that the sign of r0 is unknown as well as the higher order Fourier terms of the field d f 0 dθ .
In practice it mean that the two solutions corresponding to the two signs of d f 0 dθ have to be explored.The circular source solution is insensitive to the sign of d f 0 dθ but this is not the case for a non-circular source.As a consequence the χ 2 will have to be evaluated for each sign of d f 0 dθ and the minimal χ 2 will point point to the proper sign and break the degeneracy.
3.4.1 Reconstructing the lens and source by minimizing the χ 2 .

The parameters of the lens reconstruction.
The parameters that will have to be refined are the following, the Einstein radius R C , the center of the unperturbed circular potential r c , the pseudo impact parameter ¯r0 , and the higher order Fourier terms (N > 1) of the perturbative fields.Note that the first order Fourier terms of the perturbative fields are conjugated (see Eq. ( 21) here for δ 0 = 0), and as a consequence there are only two independent terms.The first order terms for f 1 are directly inferred from the first order terms of d f 0 dθ by conjugation.

3.4.3
The source reconstruction procedure.
The source is represented by the linear combination of 2D Gaussian polynomials (see Alard & Lupton (1998)).The 2D amplitude in the source plane represented by the coordinates (x s , y S ) is, The modelisation of complex sources which are frequently encountered in gravitational lenses requires a higher order for the expansion in Eq. ( 25).In practice to avoid the source-lens degeneracy the initial source expansion should be of low order, but should be increased to reach a final higher order with typically N = 8 or N = 10.An initial estimation of r 0 is obtained by estimating the residual to the fit of a circle to the arc (see Sec 3.3.1).Once a first fit of the lens has been performed r 0 is evaluated directly by computing the size of the actual source solution.This process of refining r 0 is iterated for each steps of the fitting procedure.

3.4.4
Fitting the image and estimating the χ 2 .
For a given model of the perturbative fields to Fourier order the use of the lens equation allows (Eq.( 3)) allows to transport each component of the source model (Eq.25) to the lens plane.Each corresponding image of the source component is then convolved with a Gaussian PSF model.The set of convolved images is then fitted to the arcs system images using the linear least-square method.This fitting procedure allows the evaluation of the χ 2 corresponding to the specific model for the perturbative fields.
Starting from the initial circular source guess we will make a search for the minimal χ 2 solution.The minimization of the χ 2 in the parameter space of the lens (see Sec. (3.4.2) ) is performed using the simplex method (Nelder & Mead (2006)).For each set of paramaters the χ 2 is evaluated according to the procedure described in Sec.(3.4.4).Due to the degeneracy on the sign of d f 0 dθ the search for the minimal χ 2 is conducted for both negative and positive signs.The choice between the 2 solutions with different signs is made by selecting the solution with minimal χ 2 .

Avoiding the source lens degeneracy problem.
The source lens degeneracy may be present due to the interaction and degeneracy between higher order source and lens terms.As shown in Alard (2017) a solution to this issue is to consider the minimal solutions.These minimal solutions are among the class minimal χ 2 solutions those with the lowest possible order in the Fourier expansion of the fields.It is demonstrated in Alard (2017) that these minimal solutions can be constructed by starting the search for the solution to the lowest possible order.A typical number to reconstruct the guess solution is N = 2 or N = 3 for the Fourier expansion of the fields.This will be also the order at which the solution is reconstructed for the first step of the refinement.Once the sign degeneracy is broken the order of the solution is increased progressively to reach the highest order for the solution (typically N = 4 or N = 5).Note that in similarity to the reconstruction process of the guess solution the size of the simplex corresponding to the higher order terms is of decreasing amplitude.Here the same rule applies, with a decrease in the amplitude proportional to order.It is possible to relax this constraint for some specific case, however it is better to try using this modulation of the higher order terms first.
3.6 Application of the method to a set of numerical simulations.
In this paragraph we will demonstrate the refinement of the guess solution using a set of 4 numerical simulations.3.6.1 Description of the numerical simulation.
The four numerical simulations corresponds to a cusp like configuration for 4 different potential types.The different potential models are a second order lens potential (elliptical) a fourth order lens potential (quadratic) and an external shear.The combination of the elliptical and quadratic potential with the external shear give the four potential models (see Table 3.6.1)and the four corresponding arcs simulations A i=1,n .

The potential models in the simulations
The potential models are based on the NFW potential described in Eq. ( 28) in (Alard (2007)), and (Bartelmann (1996)).The iso-contour of this potential are slightly distorted from circularity.The distortion is expressed using Fourier series to order N.The coefficients of the Fourier series representing the distortion are represented by η i .The lowest order corresponds to an elliptical distortion, N = 2 (2 coefficients), and the highest order in this simulation corresponds to N = 4 (quadratic).
The shear is external and corresponds to the effect of mass situated outside the Einstein circle.The shear is related to order 2 terms in the potential and as a consequence to order 2 Fourier series terms in the potential expansion.Using Eq. ( 8) for an external shear of order 2 he generic expression for the fields is, For the two numerical simulation with shear we have S 0 = S 1 = 0.05.

Building the simulations.
For each of the potential models A i the image of the source is reconstructed using ray-tracing.
The source in this simulation is complex and quite distant from a round source (see Fig. ( 2) lower panel for an image of the source contours).The image of the source was convolved with a Gaussian PSF with a FWMH of 2.5 pixels.The final image was produced by adding Poisson noise to the convolved image.
3.6.4Elliptical second order potential with no shear.
This simulation corresponds to case A 1 in Table (3 ).The refinement process starts by exploring the two solutions corresponding to each sign of the field d f 0 dθ (see Sec (3.5)).In this first iteration the Fourier order of the potential is 2. From the two solutions the solution with minimal χ 2 is retained.This solution is refined by increasing the Fourier order of the potential to 3 and subsequently to 4 to avoid the source degeneracy problem (see Sec. (3.5.1)).The refinement process is stopped at order N = 4.Not that for minimality reason the reconstruction could have been stopped at order N = 2 in this case.However in order to provide a comparison of the different lens reconstructions in a self consistent scheme all lens are reconstructed to order 4. The reconstruction of the lens to a higher is also interesting since it demonstrates the stability of the solution at higher order.All solutions were built using an eight order Gaussian-polynomial solution for the source (See Sec.(??)).The distortion vector corresponding to simulation A 1 is, η = [−0.05,0, 0, 0, 0, 0] The final result is presented in Fig. 's (2) and (3).
3.6.5Elliptical second order potential with shear.
The simulation A 2 is similar to simulation A 1 except for the presence of an external shear (see Eq. ( 27)).The construction of the guess solution and of the final refined solution follows the same process that was already described for simulation A 1 .The final result of the refinement process is presented in Fig. 's (4) and (5).3.6.6Fourth order potential with no shear.
Here simulation A 3 corresponds to an NFW potential with a fourth order Fourier distortion of the iso-contours .The corresponding distortion vector η i i define in Eq. ( 26).Here again the construction of the guess solution and of the final refined solution follows the same process that was already described for simulation A 1 .For simulation A 3 we have, η = [−0.05,0, 0.01, 0, 0, −0.005] (29) The final result of the refinement process is presented in Fig. 's ( 6) and ( 7).3.6.7 Fourth order potential with shear.
Simulation A 4 corresponds to simulation A 3 with the addition of an external shear(see Eq.   5), ( 7) and ( 9) that all configurations and even the more complex were properly reconstructed using the automated method presented in this paper.No assumptions were made about the arcs themselves and the procedure produced the lens and source modelling directly from the images.All results were produced using the same universal singular perturbative approach which suggests that an automated analysis of a large data set will now become possible with this method.Interestingly the universality of the modelling will also render possible a statistical analysis of large data sets.

CONCLUSION.
In this paper systems of gravitational arcs of increasing complexity were analysed and successfully reconstructed.Even in the case of a fourth order lens potential and a strong external shear an accurate reconstruction of the potential and source was achieved.The same complex highly non symmetrical source model was used for all the simulations.This analysis demonstrates that the method developed in this paper has the ability to perform automated reconstructions of gravitational arcs systems.The reconstructions are direct due to the na- MNRAS 000, 1-?? (2023) Fig. (1) (upper left panel).The details of this numerical simulation are available in Sec.(3.6.1).The different steps of the reconstruction are illustrated in the upper right panel with the circle fitting procedure, and in the lower panel with the first guess for the perturbative fields based on the round source model.

Figure 1 .
Figure 1.Upper left panel: the simulated image corresponding to simulation A 1 .Upper right panel: the reconstructed image.Lower left panel: the first guess estimation for the field d f 0 dθ (black cross) superimposed to the true solution(red cross).lower right panel: the first guess estimation for the field f 1 (black cross) superimposed to the true solution(red cross).
.6.1).The process of building this guess solution was described in Sec.(3.3).The resulting guess solution for this simulation is illustrated in Fig. (1).Starting from this guess solution the solution is refined in a series of steps with each time an increase of the order for the Fourier series expansion of the potential.The refinement operates according to the procedure described in Sec.(3.4.1

Figure 2 .
Figure 2. Upper left panel: a comparison between the reconstructed fields (black line) and the true solution (red line) Upper right panel: The amplitude of the Fourier series expansion of the fields as a function of the order.In both figures the field d f 0 dθ is represented as a dotted line, while the f 1 is represented by a continuous line.The lower left panel represents the true source while the lower right panels represents the reconstructed source.

Figure 3 .
Figure 3. Upper left panel: the image of the arc system corresponding to simulation A 1 .Upper right panel: the reconstructed image obtained after the final iteration.Lower left panel: subtraction between the true and reconstructed image.Lower right panel: an histogram the noise normalized by the theoretical noise expectation.Note that the histogram (red line) is quite close to the Gaussian expectation (black line).

(
27)).Here again the construction of the guess solution and of the final refined solution follows the same process that was already described for simulation A 1 .The final result of the refinement process is presented in Fig.'s (8) and (9).

3. 6 . 8
Result of the simulations.It is clear form the results of the simulations see Fig's (3), (

Figure 4 .
Figure 4. Same as a Fig. (2) but this time for simulation A 2 .

Figure 7 .Figure 8 .
Figure 7. Same as a Fig. (3) but this time for simulation A 2 .

Table 1 .
The four types of simulations.