Recovery of the Two-Dimensional Ion Distribution Function in a Magnetic Mirror from Measurements of Collective Thomson Scattering Spectra

A method is proposed for tomography of the distribution function of energetic ions that are adiabatically trapped in an open magnetic trap, according to the diagnostic data by the method of collective Thomson scattering. This method is based on measurements of the scattering spectra from successive plasma cross sections corresponding to different values of the magnetic-field strength along a single line of force. It is shown that the problem of restoring the ion distribution function in the velocity space from the measurement data in this situation is reduced to an integral equation of the first kind that allows an analytical solution. Several ways to construct exact and approximate solutions of the resulting integral equation are considered.


INTRODUCTION
Scattering of high-power probing millimeter radiation by thermal fluctuations of magnetically active plasma (the so-called collective Thomson scattering) allows one to obtain information on the velocity distribution of plasma ions with good spatial and time resolutions [1]. The possibilities of analyzing scatteredsignal spectra for diagnosing the distribution function of fast ions in toroidal magnetic traps were demonstrated experimentally on the TFTR [2], JET [3], TEXTOR [4][5][6][7], ASDEX-U [6,8] tokamaks, on the W7AS stellarator for diagnosing the temperature of thermal ions as well as the lower-hybrid instability of plasma [9][10][11], on the LHD stellarator for diagnosing the distribution of both fast and thermal ions [12], and, more recently, on the latest W7X stellarator for thermal-ion temperature diagnostics [13]. Along with optical methods and neutron and gamma-ray spectroscopy, collective scattering of millimeter radiation is one of the main methods for diagnosing the distribution function of fast ions in tokamaks [14]; in particular, this method is considered as the main technique for diagnosing thermonuclear alpha particles in the ITER tokamak reactor [15].
Recent advances in the methods for confining thermonuclear plasma in existing open magnetic traps [16][17][18][19][20] and planning of new physical tasks for nextgeneration open magnetic traps [20,21] have led to the fact that methods for diagnosing plasma in open magnetic traps are now being actively improved. An important possibility for diagnosing the distribution function of energetic ions is the implementation of the method for recording the spectra of collective Thomson scattering, which was tested in toroidal magnetic traps [22], in an open magnetic trap.
Following paper [22], we note the main features of the collective Thomson scattering in plasma of a largescale open magnetic trap using the example of the largest such GDT (gas-dynamic trap) installation that currently operates at the Budker Institute of Nuclear Physics (Siberian Branch, Russian Academy of Sciences). In this setup, plasma is divided into two fractions: background plasma that is confined in the gasdynamic regime and high-energy ions resulting from the bombardment of the background plasma by beams of neutral atoms with energies of 20-30 keV [23]. Energetic ions, first, make up a significant part of all plasma ions (20% in the central cross section and up to 100% at the points of turn); second, they are held in the trap in a collisionless adiabatic mode. The frequency of collisions of fast ions with electrons and ions of background plasma is much lower than the frequency of bounce oscillations between the points of turn in the magnetic field. In contrast to toroidal systems, the magnetic field modulus in an open trap changes greatly (e.g., in the GDT installation, the mir-

PLASMA DIAGNOSTICS
ror ratio exceeds 30); therefore, the distribution function of energetic ions is substantially different in different cross sections of the trap [22]. This, on the one hand, imposes additional restrictions on the spatial resolution of the diagnostic method and, on the other hand, can be used to restore the particle distribution function.
Let us recall that the plasma diagnostics by the collective scattering consists in the recording of the spectral power density of scattered radiation , which in turn is proportional to the spectral density of plasma density fluctuations , where ω = and are the differences of the frequencies and wave vectors of probing and scattered radiations. In the experiment, the k direction is usually fixed by the geometry of scattering, while ω changes (it falls within the spectrum-analyzer band). Therefore, we have a one-dimensional scattering spectrum that directly allows restoring of the one-dimensional ion distribution function in the velocity projections to the k direction: where is the local function of the ion or electron distribution over the three-dimensional velocity space, which is calculated in the scattering volume. For the experiment on the GDT installation that was proposed in [22], the scattering angle and the band of analysis with respect to ω are chosen so that the spectral density of fluctuations is determined by the onedimensional distribution function of fast ions.
Hereafter, we assume that for a particle that rotates in a magnetic field, all directions of the velocity across the magnetic field are equally probable; therefore, an ensemble of such particles can be characterized by a two-dimensional distribution function in the transverse and longitudinal velocity-vector components with respect to the magnetic field. Taking advantage of this natural symmetry, we can rewrite the expression for the one-dimensional distribution function as (1) where θ is the angle between the k vector and the magnetic field direction. Measuring the scattered-signal spectrum at a given scattering angle does not allow one to restore the full distribution function , but such a restoration is possible in principle using the two-dimensional function , i.e., if we perform measurements for a certain set of scattering angles. This approach to the recovery of the distribution function is close to the classical tomography problem, which can be solved in this case via standard mathematical methods for solving ill-posed inverse prob- lems. This method was proposed and theoretically analyzed in [24], demonstrated on synthetic data [25], and successfully implemented when processing the actual experimental data [14].
However, a wide variation of the scattering angle is not always possible and requires a rather sophisticated geometry of the lines of injection and reception of microwave radiation. In the case where measurements can be performed at points that correspond to different values of the magnetic field along the same line of force, another approach is possible: we can take advantage of the fact that if the electron distribution function is recalculated between these points in a known way, we can also restore the full distribution function by comparing the scattering spectra at these points. In this paper, we mathematically formulate the problem of such a recovery, show that such a problem can be reduced to finding solutions to an integral equation of the first kind, and demonstrate ways to solve this equation.

THE INTEGRAL EQUATION
Let us consider formula (1) for an arbitrary point along a magnetic-field line of force: where is the distribution function in the cross section where the magnetic field strength is , and the local mirror ratio R plays the role of the coordinate along the magnetic field line of force. Hereinafter, we use the notation for the velocity in a specified cross section R and for the velocity in the central cross section , which corresponds to the minimum magnetic field value. For an ion that performs collisionless motion in a smoothly inhomogeneous magnetic field, two integrals of motion can be introduced The first integral corresponds to the conservation of the magnetic moment, while the second corresponds to the conservation of energy with allowance for the electrostatic ambipolar potential . The characteristic values of the ambipolar potential drop along the line of force are determined by the target-plasma temperature [18,26]; therefore, it can be ignored when considering the dynamics of energetic ions. As a result, the relationship between the velocities and for fast ions can be found in the explicit form (2) By virtue of the Liouville theorem, the phase volume is conserved; therefore, the ion distribution function in the given cross section R is expressed through the distribution function in the central cross section of the trap in the form: where the velocities at the minimum of the magnetic field are expressed through the velocities in the cross section R according to transformations (2).
For simplicity, we choose the longitudinal scattering geometry and take the integral over . We obtain Further, we change to the velocities at the magneticfield minimum in the integrand: Formally, the integration is performed over the domain of , into which the half-plane of , transforms. However, the δ-function is equal to zero everywhere beyond this domain; therefore, we can extend the integral to the entire space of velocities . In addition, owing to the symmetry of the distribution function for trapped ions that is expected in an open trap, we can consider only the domain of for ; in this case, the subintegral δ-function is equal to zero in the domain of . Taking this into account, we can rewrite the integral in a simpler form: Taking the integral over , we obtain , , This is the main relationship that will be analyzed under the assumption that we know the continuous function from the results of the plasma diagnostics. Of course, when using the data from an actual experiment to construct such a function, we will have to use an interpolation over a discrete set of measurements.
For the compactness of the subsequent calculations, we introduce a new notation where is any typical velocity introduced only for the normalization. As a result, we come to the integral equation (4) The mathematical statement of the recovery problem consists in finding the function using a specified function . In this case, we consider that the variables change in the interval . The function is a real positively defined piecewise-smooth function, the integrals of which over any finite domain are limited. It should be noted that the real distribution function must also vanish starting with certain values of the variables because of relativistic limitations. However, for the convenience of calculating analytic integrals, we will not require the fulfillment of this condition, but instead reserve the freedom to impose restrictions on the law of decay of the function for .
3. THE FOURIER TRANSFORM Since all the arguments in Eq. (4) are given on the positive half-axis, it is natural to use the Fourier cosine transform to find a solution. We define (5) A necessary condition for these equalities is the convergence of the integrals for all values of . After transformation (5), Eq. (4) is rewritten in the form As is seen, the Fourier cosine transform preserves the integral character of the equation, but reduces it to a normal form that admits a standard solution. This solution is presented in the next section; however, before proceeding to it, let us make a digression of a methodological character.
The functions and are fundamentally defined only at positive values of their parameters. Let us forget for a while about this and extend all our dependences to the domain of ; however, it is considered as before that . We then can apply a complex Fourier transform in the variables to Eq. (4): In this case, Eq. (4) takes the form Next, we take into account that the integral is taken only for positive y values and come to For , this relationship gives a simple algebraic connection between the Fourier transforms of the preset function G and the desired function f. However, for , we obtained an additional restriction on the function , which must be taken into account when extending the definition of its original to the domain of negative arguments. For the essence of this restriction to be understood, we can transfer the uncertainty and restrictions of the function G at to the uncertainty and restrictions of the function in the domain of . Let us take advan- tage of the fact that and are real, i.e., and , where the asterisk designates the complex-conjugation operation. Without a loss of generality, we consider that the distribution function is redefined symmetrically in the nonphysical domain: . The complex Fourier transform then turns into cosine transform (6), and . Using these conditions, we have This equation is equivalent to Eq. (6). Thus, if the function is defined correctly in the domain of , then, using a complex Fourier transform, we can easily solve Eq. (4). However, the problem of correctly redefining is anyway reduced to the solution of Eq. (6), which is generated by the Fourier cosine transform.

SOLUTION VIA THE MELLIN INTEGRAL TRANSFORM
Equation (6) is formulated for an unknown function of two variables. Despite its simple form, we were unable to find this equation in mathematical reference books. Nevertheless, its solution can be constructed based on the technique that was developed for solving one-dimensional linear integral equations of the first kind [27].
Let us introduce the Mellin integral transform [28] The integral in the second expression is taken in the sense of the principal value, while the parameter σ belongs to the interval , whose boundaries are determined by the convergence of the integrals The convergence of these integrals for some at is the criterion for the applicability of the Mellin transform for the function .
Let us multiply the left and right sides of Eq. (6) by , integrate over ξ and w from zero to infinity, and obtain (7) The left side of Eq. (7) contains the two-dimensional Mellin transform for the function : the inverse transform to which has the form The range of changes of the parameters , will be determined after the finish of the transformation of the right side of Eq. (7). The expression in square brackets in (7) is the Mellin transformation for the kernel of integral transform (6). We transform this expression by moving to the new integration variables: The expression in square brackets can now be calculated using the tabulated Mellin transforms for trigonometric functions [28]: (9) where denotes the gamma function, and the calculations are correct for . After the above transformations, expression (7) takes the form where (11) Let us determine the required restrictions on , . We depart from the existing constraints for integral transform (9): Re(s ξ ), (0, 1). For Mellin transform (8) for the function , this range of the parameters requires that the positively defined function have no singularities at and and decay at infinity no slowly than . The conditions for Mellin transform (11) for the function are determined somewhat more difficultly because it is not positively defined. Let us express the desired integral through : The integral over is taken explicitly and corresponds to the Mellin transform for [28]. The convergence of the integral with respect to v requires the additional condition ; if it is met, then taking the conditions of into account, the integral over x converges, provided that has no singularities at and decreases more rapidly than for . When considering the analogous transform with respect to y, we obtain that the convergence requires that have no singularities at and decay more rapidly than for ; in this case, additional limitations on do not arise. We note that the conditions imposed on the distribution function are fulfilled automatically, if the natural physical requirement of finiteness of the number of particles and the total kinetic energy in the system is imposed. Next, it can be shown that if the above conditions are met for f, the analogous conditions are also met for the function G.
Finally, for the applicability of the Mellin transforms, it is sufficient to require that the constants in the Mellin transforms belong to the following intervals: The solution of Eq. (6) can be represented in the form of the following sequence of linear transformations: The numbers of the formulas that describe transforms are indicated above arrows. This chain is expressed in the explicit form as Changing the order of integration allows one to obtain the solution in the standard form for a linear operator in the form of an integral convolution (14) with the kernel The latter expression can be simplified, since, as already noted, the integral over is taken explicitly and corresponds to the Mellin transform for ; the resulting trigonometric functions are cancelled, and the following expression remains: (15) This expression can be considered as the definition of a function of two but not four variables, since depends only on the combinations of and . In the future, we will not need this self-similarity property, however, it allows tabulation of the kernel using a two-dimensional table or interpolate it using a function of two variables. This can be the basis of an efficient practical way to calculate the twodimensional convolution (14).
Let us prove that the resulting inverse integral operator actually allows restoring of the function . Operator (14) with kernel (15) is applied to the function of the most general form (4): Next, we change the variables, , and change the order of integration so that the integrals over and are taken last of all, In the resulting expression, the integral over w is the known tabulated integral [28]: (17) that converges at , , . Substituting (17) into the initial expression (16), we obtain that all gamma functions are cancelled and from the remaining part, the direct and inverse Mellin transforms are taken [28], which form a delta function: Thus, we have shown that the recovery of the original function is correct within the limits of the applicability , .
y f x y ds ds dx dy f x y x x y y dx dy f x y v of the Mellin transforms. In Appendix, one can additionally find several particular examples of the restoration of the distribution function, which appeared on the way to finding a general solution and seemed to be quite interesting.
It should be noted that the mathematical proof of the existence and uniqueness of solutions to Eq. (6) for arbitrary G is quite laborious, since it is an equation with an alternating kernel [27]. However, the fact that as G, we take only the functions that are transforms for some quite rapidly decreasing function f, which arises from the physical problem, is sufficient for our solution to exist and be determined by relationships (14) and (15).

ORTHOGONAL-POLYNOMIAL EXPANSION
OF THE SOLUTION Although we managed to present an analytical expression for the inverse operator in quadratures, nevertheless, this solution is largely formal and inconvenient for an actual application. Therefore, here, we propose another way to invert Eq. (4).
Using the given function , we can restore the discrete set of functions defined as It is easy to show that (18) i.e., these functions describe the dependences of the moments of the distribution function in the variable y on the variable x (recall that x and y are, respectively, the kinetic energies of the longitudinal motion and cyclotron rotation of a particle). Linear combinations of the functions can be composed so as to provide the formation of some polynomials that are orthogonal on the interval. Let us consider Laguerre polynomials as the basic example [29]: with the orthogonality condition Here, denotes a binomial coefficient. We introduce the notation for the formal Laguerre polynomial , in which the replacement is performed: . In view of (18), we obtain a set of integral equations for (19) Here, the variable x is merely a parameter, and the variable y can be used to expand the distribution function in terms of the full system of Laguerre polynomials: (20) where are the expansion coefficients that are determined independently for each x value. Owing to the orthogonality of the Laguerre polynomials, Eqs. (19) assume the form Thus, we have constructed the solution of Eq. (6) in the form (21) The resulting series has an important property that ensures its effectiveness in solving the physical problem. The partial sum of the first N terms over n gives the function that has the same set of the first N moments in the variable y as the desired function has; thus, if we substitute into relationships (18), they then are valid strictly for for all x. Correspondingly, the remainder term has the zero first N moments. For example, considering only the first three terms yields the correct mean value and variance of the distribution in y, which depend on x as a parameter.
Let us discuss the convergence of the resulting series. To do this, it is sufficient to consider the analogue of Parseval's theorem for expansion (20), which follows from the conditions for the orthogonality of Laguerre polynomials, . To estimate the rate of convergence, let us consider a function in the form (22) The coefficients of expansion (20) for this function are It can be seen that when the convergence condition is met, the amplitude of the coefficients decreases exponentially with increasing n. That is, for an infinitely smooth function, the coefficients of series (21) decrease exponentially. In other words, if the series converges, it converges quickly. In the case where there is a discontinuity in the value of a function or its first derivative, the rate of convergence can also be estimated using particular examples. Let us consider the expansion in terms of Laguerre polynomials for the model functions and that have discontinuities of a value of the function and the first derivative, respectively; here, Figure 1 shows the partial sums of the first 11 and 101 terms of the series for the function (on the left) and (on the right). It can be seen that for a function with a discontinuity in its value, taking even one hundred terms of the series into account gives no satisfactory convergence; however, for a function with a derivative discontinuity, a satisfactory convergence is achieved ( )   Note that the divergence of series (21) for distribution functions with power-law tails follows from the scheme of constructing our solution itself, which requires the convergence of the set of integrals (18). We cannot use terms of the series with for restoring the distribution functions with power-law tails . Nevertheless, even in cases where our series does not converge formally due to the fact that high moments of the original function tend to infinity, we can use a finite number of finite terms of series (21) for a satisfactory approximate solution of the original problem. From the physical point of view, this is not surprising, since all the lower moments for the approximate distribution function are guaranteed to coincide with exact values for the restored function f. The procedure for constructing an approximate solution using a finite number of series terms is discussed in the next section.

RESTORING THE FUNCTION FROM A FINITE NUMBER OF TERMS OF THE SERIES
In order to clarify the idea developed below, let us start with a simple special case. We have certain arbitrariness when choosing a system of orthogonal polynomials. In particular, we can use polynomials , where does not depend on y, instead of . In view of a rescaled orthogonality relationship, it is easy to obtain a set of various equivalent expressions for series (21) with an additional free parameter β: This series inherits the main property of expansion (21)-the partial sum of the first N terms gives an approximate function that exactly reproduces the first N moments with respect to the variable y, that is, in this case, the remainder term includes the zero first N moments.
It can be seen that by selecting β, we can provide the convergence of the resulting series for any exponentially decreasing function f, including those for which the initial series (21) diverges. Moreover, it is not difficult to see that all functions of the form (22) will be exactly described by one term of series (23) if . In a real problem, the function f is unknown a priori; moreover, the asymptotic behavior of this function is most likely not exponential for . Nevertheless, the freedom to choose the parameter also allows optimizing the convergence of infinite series (23) for f and decreasing the error of the approximate solution , which is related to the rejection of the remainder term. Such a rejection in an actual problem is primarily due to the natural limitations of the experimental data, which cannot be infinitely differentiated with the hope to extract some information from them. This raises the problem of optimizing the restoring of the distribution function in a situation where only a finite number of series terms are available. We can propose the following approach to its solution. The first N terms of the expansion are taken into account, and by choosing β, we ensure that the th term of the expansion is equal to zero. This implies solving the algebraic equation (24) for β at each x value. The obtained approximate function will consists of N expansion terms and precisely reproduce the first moments of the distribution function with respect to y. Bearing in mind the meaning of the variables x and y, Eq. (24) can be interpreted from the physical point of view as the equation for the effective transverse temperature, which most accurately describes the distribution of particles with a specified longitudinal energy.
The described method can be improved if generalized Laguerre polynomials that are defined as where denote the generalized binomial coefficients [30], are taken as the system of orthogonal poly-   (25) where is considered as the expression obtained from the polynomial as a result of the formal substitution . The obtained series converges in a wide range of the parameters α and β for the function that decays rather quickly at infinity. This allows us to optimize the approximation of the distribution function by a finite number of terms of the series by varying these parameters; these parameters can be varied independently for each x value.
As before, the partial sum of the first N terms of series (25) exactly reproduces the first N moments of the desired function at any α and β. Let us demand that the th and th terms of the series be strictly equal to zero. This will lead to the algebraic system (26) which determines the parameters and independently for each x value. The approximate function that corresponds to this choice of parameters will exactly reproduce already moments of the distribution function.
The expansion in terms of generalized Laguerre polynomials allows one to take the specifics of the particle distribution in an adiabatic magnetic trap into account, namely, the distributions with an empty loss cone. Recall that y corresponds to the kinetic energy of the transverse motion (rotation) of a particle. Therefore, distributions with a loss cone are well modeled by functions of the form that can be described exactly by one term of series (25). It can be expected that when restoring the distribution function with a depleted ion distribution in the loss region from actual data, expansion (25) will provide the preset accuracy with fewer terms of the series than expansion (23), which relies on "less physical" basis functions.  Let us consider the examples of how the recovery of the distribution function works through the expansion in terms of Laguerre polynomials with one varying scale β and through generalized Laguerre polynomials with two varying parameters α and β. The verification procedure is the same all the time: using the analytically defined function , we calculate the diagnostic response , calculate the moments , and then apply solutions (23) or (25) with a finite number of terms of the series. We always consider three terms taking the forcedly zeroed ones into account; i.e., the restored approximate function contains at most two terms if one expansion parameter is varied, and only one term if two expansion parameters are varied. Even under such strict conditions, Eqs. (24) or (26) guarantee that the first three moments, i.e., the number of particles, the total energy of cyclotron rotation, and its variance, are restored exactly for each group of particles having a given longitudinal velocity. Figure 2 shows the results of restoring the distributions of the form . In this case, the multiplier that is responsible for the dependence on x is restored exactly by the finite sum of the series; in this case, the parameters α and β are the same for all x. Figure 2a refers to the elementary case of an anisotropic Maxwell distribution, for which the first term of expansion (25) gives an exact answer for and (the same graphs on the left and right). If we fix the "wrong" transverse temperature and optimize the parameter α, the first two terms of the expansion then provide a relative accuracy at a level of 5% in the entire "thermal region" (the graph at the center, the relative accuracy is determined by the maximum relative shift of the level lines). The case (b) corresponds to a Maxwellian distribution with a loss cone, for which the first term of expansion (25) gives an exact answer for and (the graph on the right). The other two methods (on the left, fitting the transverse temperature β at and at the center, fitting α at the forced incorrect temperature β = 1) allow satisfactory recovery of the distribution function over only two expansion terms. The case (c) also corresponds to the Maxwellian distribution with a loss cone, but it has a more complex form. Now we already cannot accurately reconstruct the distribution function by a single term of the series. It can be seen that the optimization simultaneously in two parameters, α and β, which leads to a simple formula ≈ (graph on the right), gives a better approximation than the two-term approximations that are obtained as a result of optimizations in each parameters separately (on the left and at the center). Thus, using the above "training" examples, we became convinced that expansion (25) with an additional free parameter can significantly improve the x . y accuracy of the recovery compared to one-parameter expansion (23).
Let us now consider an example that simulates a physical problem. Figure 3 shows the level lines of the distribution function that describes hot ions generated as a result of both stripping a neutral beam, which is injected into plasma at an angle of 45° to the magnetic field, and the subsequent Coulomb relaxation, which is the case that meets the experimental conditions on the GDT installation. The initial distribution function is specified either by an anisotropic Maxwell distribution that is shifted to the point corresponding to the velocity of neutral atoms (the graph on the left), (27) or a similar distribution that is truncated on the side of high energies (the graph on the right) (28) The second distribution simulates the features of relaxation of the distribution function due to Coulomb collisions with background plasma particles, in which the diffusion in the velocity directions and energy loss due to friction are the dominant effects. Unlike the effects discussed above, these distributions are no longer factorized; therefore, when restoring them, we performed optimization with respect to α and β separately for each x value. The results of the recovery are shown in the graphs with dashed level lines. The restoring accuracy of the function using only one term of series (25) is approximately 1%. It can be seen that for a function with a discontinuity, the recovery is somewhat worse, but its accuracy is still quite high. As an example, the most difficult parameter to recover by this method is the position of the discontinuity that corresponds to the maximum energy ; it is restored with an accuracy of ~5% when determining the maximum energy at a level of the distribution function of 0.01. We also note that the approximation obtained from a small number of terms of the optimized series for a function with a discontinuity is positively defined, in contrast to the result obtained by summing a large number of terms of an unoptimized series (Fig. 1). We recall that the method automatically provides correct values for the number of particles, transverse energy, and its variance at each value of the longitudinal particle velocity.
Hence, by summarizing our experience in applying this method to various distribution functions, we can formulate the following empirical conclusion. If the choice of α and β ensures that two consecutive terms of series (25) are equal to zero, then the sum of all subsequent unaccounted terms of the series will be small ( ) in comparison to , or, in other words, the sum of the first nonzero terms will repeat well the desired function. This statement is also true for individual func-N f tions with "pathological" features, e.g., for functions with power tails, for which all the higher moments, beginning with a certain one, diverge. 7. CONCLUSIONS In this paper, we have studied the possibility of reconstructing a full two-dimensional velocity distribution function of energetic ions in an open magnetic trap from measurements of collective Thomson scattering spectra in the same scattering geometry, but in different trap cross sections that correspond to different magnetic field values. It is seen that under the assumption of the adiabatic nature of the ion motion in the trap, such a series of measurements is sufficient to completely restore the distribution function.
The problem of reconstruction is reduced to the solution of an integral equation, for which two methods of solution are proposed. Both presented analytical methods, through the Mellin transform and through an expansion in terms of orthogonal polynomials, allow realization in practice. However, the method of expansion in terms of orthogonal polynomials appears somewhat more preferable from the standpoint of the physical validity and ease of implementation. When using the expansion in terms of orthogonal polynomials, each term of the series is responsible for restoring the corresponding moment of the distribution function. Naturally, the search for a solution in the form of a series also has its drawbacks, which include the need to calculate derivatives of increasingly higher order of experimental data for calculating terms of the series. In this case, the problem arises of optimizing the procedure for restoring the distribution function when only a finite number of terms of an infinite series are available from the experimental data in principle. The possibilities of such optimization are considered in this paper on the example of a case where only three terms of a series are available, which requires calculating the first and second derivatives of the experimental data. The possibility of a very good recovery of model distribu-tion functions due to the choice of the scale and order of generalized Laguerre polynomials is shown.

APPENDIX
Here, we consider examples of the reconstruction of the distribution function via the Mellin transform, which allows a fully analytical calculation. These particular examples appeared during the verification of formulas (14) and (15) proposed in the article before we obtained the general formal proof presented here. They can also be useful for tests during the practical implementation of the proposed method.
I. Delta function. Let us consider the distribution function that corresponds to a monovelocity ion beam, The function assumes the form Let us apply inverse (14) to the function . This results in the following chain of expressions:   Fig. 3. The recovery of the ion beam distribution function using the first term of series (25) for the optimization in both parameters, α and β. The left and right graphs correspond to smooth distribution (27) and distribution (28)   The integral over ξ corresponds to the already used tabulated integral (17) and is taken under the condition of and ; the latter equality follows from the Mellin transform for the delta function (the rigorous proof of the fact that the obtained integral is a delta function can be found in [31]).
II. Anisotropic Maxwellian distribution function. Let us consider the distribution function The function that corresponds to it has the form Further, the tabulated integrals [28] are used: (A.1) (these formulas are valid for and ). As a result of applying inverse integral operator (14) to the function , we obtain During calculations, the relationships and were used as well; the latter equality follows from the known Mellin transform for the exponent [28] which is valid for ; or . As is seen, the latter restrictions lie in domain (13).

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Theimages or other third party material in this article are included in thearticle's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.