Iterative algorithm for interferometric retrieval of plasma density in case of considerably inhomogeneous objects

Interferometry is a widely used active diagnostic method for measurements of physical properties of plasmas. In case of axial symmetry of probed objects an Abel inversion can be used to retrieve plasma density profiles from phase shifts reconstructed from interferograms. This approach is based on the assumption of the diagnostic beam propagating along a straight line. However, it is well known that in case of inhomogeneous media the refraction process affects the diagnostic beam trajectory resulting in an inaccuracy of the retrieved density structure. In order to deal with this unfavourable effect a more sophisticated approach needs to be employed. In this paper a special iterative algorithm is proposed to deal with this issue. This algorithm turns the inversion procedure into series of iterations, where the appropriate distribution of plasma density is found by following the diagnostic beam actual trajectory during its propagation through the inhomogeneous medium. The assumption of axially symmetric plasma distribution still applies. The respective trajectories are calculated using the ray tracing method. This method also allows to use the paraxial wave equation. This iterative algorithm was tested on simulated data with different configurations of plasma density proving its functionality. Results from the simulated data analysis show that using this approach the effect of refraction can be fully compensated and plasma density is thus recovered accurately. Comparison with the results obtained only by the Abel inversion as such is provided for illustration.


Introduction
Optical interferometry belongs to important active diagnostics for laser produced plasma studies. Its main application is for retrieval of plasma density spatial distribution from the phase shift reconstructed from interferograms. In experiments involving the laser produced plasma or other optically transparent objects in general, certain symmetries can be encountered. In this 14th Kudowa Summer School "Towards Fusion Energy" IOP Conf. Series: Journal of Physics: Conf. Series 1197 (2019) 012002 IOP Publishing doi:10.1088/1742-6596/1197/1/012002 2 paper, axially symmetric objects are considered. The widely used approach towards retrieval of plasma density from the phase shift obtained by interferometry of axially symmetric objects is based on the Abel inversion technique [1]. This technique assumes the propagation of the probing beam along a straight line through the measured object. In practice, however, there are some important drawbacks of this assumption which should be taken into account. As the probing beam is propagating through the plasma, which is in majority of practical cases an inhomogeneous object, its wavefront gets distorted as well as the distribution of its intensity across the cross section. These features are connected with the refraction effect and dealing with their consequences will be a main goal of this paper.
At this place it should be noted that a generalization of the classical optical interferometry called Complex Interferometry was developed in the past to retrieve from just one complex interferogram up to three quantities generated by the optically transparent object under investigation: the probing beam phase shift, the change of the probing beam amplitude as well as the fringe contrast [2]. Some aspects of this more general technique will be employed.
Over the years several approaches aimed at resolving the effect of refraction on the optical probing were considered. In the case of axially symmetric objects computational inversion schemes were proposed [3]. In this paper, a novel approach towards retrieval of plasma density from interferograms of axially symmetric objects affected by the effects of refraction is presented. This approach was verified on simulated plasma density distributions where the effect of refraction was present. This method is based on iterative algorithm, using successive corrections to the initial estimate of the plasma density distribution obtained by the Abel inversion. The verification using simulated profiles shows a quick convergence of the algorithm and demonstrates its ability to compensate for the effect of refraction that is not accounted for using the Abel inversion alone.

Theoretical part
In order to evaluate the effects of refraction a description of the light propagation through inhomogeneous media is required. The geometrical optics approximation can be used for transparent inhomogeneous objects under the assumption that the characteristic length of changes in material properties (index of refraction) is significantly larger than the wavelength of the probing beam.
The method used to determine the light propagation through the inhomogeneous media is based on the solution of the ray equation (two possible forms) [4] The rays are expressed as dependent on parametric coordinates: τ corresponds to the ray parametrization, while s corresponds to physical length. The advantage of the ray equation using τ parameter is its correspondence to the equation of motion of a particle in a potential field. This allows to use a variety of differential schemes to numerically solve this equation. The ray equation is, in general, obtained from the solution of the Helmholtz equation using the power series in the inverse wave number (2). The assumed solution consists of the eikonal ψ and the amplitudes A m , where the zeroth element gives the geometrical optics approximation and the higher ones take into account the effect of diffraction. The second element of the power series A 1 can be used to verify the accuracy of the geometrical optics approximation using following formula γ = |A 1 |/|k 0 A 0 | which should be significantly lower than one in order for geometrical optics approximation to apply. It should be noted that the conditions for geometrical optics approximations put constraints on the spatial variations of the medium through which the rays propagate.
The optical path length is expressed by the eikonal which is given by integrating the refractive index along each ray that was obtained using the solution of the equation (1). The formula to calculate the optical path length is derived on page 11 in the book [4].
From the optical path length the absolute value of the phase can be obtained provided the wavelength is known. Solution of the ray equation (1) gives the trajectories of the rays propagating through the inhomogeneous medium. Using the ray positions (initial and final) the corresponding amplitude can be calculated from the conservation of energy approach illustrated in the Figure 1 using the following formula Figure 1. Conservation of energy along the ray tube dependent on the amplitudes A 0 , cross sections of the ray tube da and the refractive index n. Superscript (0) denotes the initial values.
From the phase shift given by the eikonal and the amplitude calculated from the geometrical positions of the rays the resulting scalar wave function u can be calculated.
In order to retrieve the spatial distribution of axially symmetric objects the Abel inversion is employed. This method assumes that there is no refraction and the rays along which the beam propagates are straight lines. This assumption is sufficient for objects with reasonably small refractive index gradients where practically no effect of refraction is observed. However in practice, refraction can be shown to have a detrimental effect on the accuracy of the reconstructed data when only the Abel inversion is employed. This is illustrated in Figure 2 where a propagation through simulated two dimensional axially symmetric object was used to provide the corresponding phase shift which was subsequently inverted using the Abel inversion. As the non-linearity of this task does not allow for an algebraic formula for the inversion process to be performed, computational methods needs to be employed [5]. Due to this reason an iterative algorithm based on successive approximations is used instead. For the forward propagation the ray tracing method can be used. Original density Fourier series Onion peeling Figure 2. Inaccuracy of the Abel inversion applied to a considerably refracting object.

Implementation
The sequence of the iterative algorithm steps is outlined in the list below. Typical input for this algorithm is the measured phase shift ϕ M reconstructed from an experimentally obtained interferogram. In this paper, however, such phase shift is obtained from the ray tracing simulation of the probing beam propagation through a defined object. In every iteration step all physical quantities are denoted by a superscript indicating its number. This applies to plasma density n (i) , phase shift ϕ (i) , and correction function f (0) .
(i) First estimate of the plasma density n (0) is obtained from the measured phase shift ϕ M using the Abel inversion and the counter is set to i = 0. (ii) Using the ray tracing the phase shift ϕ (i) and the corresponding amplitude are calculated. (iii) Correction function f (i) = ϕ M /ϕ (i) is obtained as a division between the measured and the calculated phase shifts. This correction function is further modified to account for the geometric distortion of the curved rays. For this, the points where the correction function is sampled are recalculated. The correction function is then recalculated to the points, where the plasma density is defined. (iv) New plasma density n (i+1) = n (i) f (i) is calculated as a more accurate approximation. (v) If the resulting plasma density tends to converge towards some fixed value fulfilling the convergence criterium ||n (i+1) − n (i) || < δ the algorithm is stopped, otherwise a jump to the step (ii) is performed.
The iterative algorithm for the plasma density reconstruction outlined above was implemented in the GNU Octave programming environment. Inputs for the algorithm are the measured phase shift retrieved from the experimentally obtained interferogram, the dimensions of the input image and the wavelength of the probing beam.
To fully compensate for the effect of refraction, the effect of geometrical distortion that results from the propagation through an inhomogeneous medium needs to be taken into account. The iterative algorithm uses the correction function f (i) which modifies the retrieved plasma density. The effect of geometrical distortion can be illustrated in the figure 3, where the rays indicate the geometrical mapping from the points where plasma density is sampled to the corresponding positions on the final wave surface, where the phase shift is recorded.
The calculated plasma density after every iteration step is numerically filtered using the convolution based smoothing. Simulation of the probing beam propagation through the measured object is outlined in Figure 3. The spatial distribution of plasma density of the axially symmetric object is given by the function that is calculated in each iteration. Final and initial wave surfaces only denote the numerical solution of propagation. In an experimentally obtained interferogram, the retrieved phase shift is determined by the propagation of probing beam through the measured object up to the detector. In these simulations the phase shift is directly determined after the beam propagates through the object.   The formula used for the refractive index of the electron plasma density is given below, where the plasma frequency is denoted as ω p .
In these formulae the electron plasma density is denoted as n e and the angular frequency of the probing beam as ω, 0 is permitivity of vacuum and m e is the rest mass of an electron.

Results
The iterative algorithm was tested on computer generated data with spatial features similar to experimentally observed plasma density profiles. The computer generated data were chosen for a possibility of comparison with the correct plasma density. At first, the 2D version of the iterative algorithm was tested. The results are shown in Figure 4. The algorithm converges in 10 iterations, depending on the set condition for conversion. Similar results of convergence were obtained for a spatial distribution with a central dip. The first series of tests can determine possible use of the iterative algorithm to compensate for effect of refraction in experimentally obtained interferograms.
The wavelength used for this simulation was 800 nm, size of the object was 200 µm. The plasma density was 1.35 x 10 26 m −3 at the top of the spatial distribution which featured a cosine shape. Under these conditions, the geometrical optics approximation is sufficiently accurate.
Propagation through a 3D object was also tested. The iterative algorithm was observed to also converge towards the original plasma electron density values. The object had a simple geometry, a ball of plasma with index of refraction below unity. Yet, these parameters had a considerable refraction effect present. The results of propagation through this 3D object with spherical symmetry in spherical coordinates are shown in Figure 5. The visible artifacts are due to remapping of the refractive index to cylindrical coordinates in the ray tracing implementation. The remapping is done using linear interpolation between sampled points. Possible way to reduce these artifacts is to increase the sampling of the both rays and spatial distribution of plasma density and then use numerical smoothing.

Discussion
It has been observed that a numerical filtration of the results of the simulation improves the convergence of the iterative algorithm. The filtration method is based on moving average smoothing and is implemented using convolution. Due to the remapping of the sampling points of the physical quantities the iterative algorithm had a tendency to numerically oscillate at the output. This has been successfully compensated for and with this improvement the algorithm was shown to converge even for slightly more complex simulated objects. The experimentally obtained data always contain certain levels of noise. The noise manifests itself as a high frequency part on top of the low frequency background. However, this high frequency part is usually filtered depending on the exact method of the interferogram data analysis. Provided the method based on the fast Fourier transform is employed the inherent low pass filtering gives smoother results automatically. Given the limits on the use of the algorithm only reasonably smooth spatial distributions of plasma density and phase shift can be considered for analysis so far.

Conclusion
In interferometric measurements of plasma density the effect of refraction on the propagation of the probing beam is known to have a detrimental effect on the accuracy of this measurement provided only the Abel inversion alone is employed. In order to deal with this feature some method which compensates the effect of refraction is required. For this purpose one iterative algorithm was proposed. This algorithm uses ray tracing to follow the actual probing beam propagation trajectory through the measured object. The zero approximation estimate of the plasma density is obtained using the Abel inversion. Subsequently an iterative process is used to find the corrected plasma density spatial distribution which fits the measured phase shift.
The iterative algorithm proposed was successfully verified on simulated axially symmetric objects in both two dimensional and three dimensional propagation of the probing beam. As this algorithm is just in its phase of initial tests and verification further development is needed. Future possibilities include improvement of the ray tracing code and iterative algorithm for different geometries.