Bayesian inference of ion velocity distribution function from laser-induced fluorescence spectra

The velocity distribution function is a statistical description that connects particle kinetics and macroscopic parameters in many-body systems. Laser-induced fluorescence (LIF) spectroscopy is utilized to measure the local velocity distribution function in spatially inhomogeneous plasmas. However, the analytic form of such a function for the system of interest is not always clear under the intricate factors in non-equilibrium states. Here, we propose a novel approach to select the valid form of the velocity distribution function based on Bayesian statistics. We formulate the Bayesian inference of ion velocity distribution function and apply it to LIF spectra locally observed at several positions in a linear magnetized plasma. We demonstrate evaluating the spatial inhomogeneity by verifying each analytic form of the local velocity distribution function. Our approach is widely applicable to experimentally establish the velocity distribution function in plasmas and fluids, including gases and liquids.

The velocity distribution function is a statistical description that connects particle kinetics and macroscopic parameters in many-body systems. Laser-induced fluorescence (LIF) spectroscopy is utilized to measure the local velocity distribution function in spatially inhomogeneous plasmas. However, the analytic form of such a function for the system of interest is not always clear under the intricate factors in non-equilibrium states. Here, we propose a novel approach to select the valid form of the velocity distribution function based on Bayesian statistics. We formulate the Bayesian inference of ion velocity distribution function and apply it to LIF spectra locally observed at several positions in a linear magnetized plasma. We demonstrate evaluating the spatial inhomogeneity by verifying each analytic form of the local velocity distribution function. Our approach is widely applicable to experimentally establish the velocity distribution function in plasmas and fluids, including gases and liquids.
The Maxwellian distribution is the most fundamental velocity distribution function, which holds in thermodynamic equilibrium states. However, it does not necessarily hold for non-equilibrium states. The ions in magnetized plasmas, which are usually far from equilibrium, are subject to other types of the velocity distribution function. The bi-Maxwellian distribution represents the anisotropy of flow and temperature in parallel and perpendicular directions to the external magnetic field 1 . The ring velocity distribution is a model of the energetic ion population that is formed immediately after neutral beam injection 2 . The velocity distribution function of α-particles has a fat tail originating from collisions of these particles with the thermal ions and electrons in the background plasma 1,3,4 .
Laser-induced fluorescence (LIF) spectroscopy is utilized to observe the velocity distribution of ions or neutral particles in low-temperature plasmas [5][6][7] . Local features in inhomogeneous plasmas, such as number density, flow velocity, and local temperature, are obtained by fitting the representation of the velocity distribution function in the absorption frequency domain to an observed spectrum. However, the explicit expression of such a function is often unclear or ambiguous. While intricate factors deform the velocity distribution in non-equilibrium states, as highlighted above, it is not always easy to tell which is dominant/negligible. Furthermore, changes in internal degrees of freedom, such as the Zeeman and Stark effects, make it more challenging to identify the explicit expression of the velocity distribution function. They deform the shape of the observed spectrum even if the velocity distribution itself is unchanged 6,8,9 . It is crucial to experimentally establish the analytic form of the velocity distribution function and its parameters without ambiguity for understanding the system of interest.
In this study, we propose a novel approach to evaluate which form of the velocity distribution function assumed for an observed spectrum is valid. Our approach is based on the Bayesian model selection [10][11][12][13] , which quantify the uncertainty in each analytic form of the function assumed for observed data. We formulate the Bayesian inference of the velocity distribution function from an observed spectrum. To demonstrate the validity www.nature.com/scientificreports/ of our approach, we apply our approach to LIF spectra locally observed at some positions in a magnetized helicon plasma. We successfully verify that each local ion velocity distribution function is Maxwellian rather than other candidates.

Results
Models. We introduce our models while briefly reviewing the relationship between the velocity distribution function and their representation in the absorption frequency domain. To observe the velocity distribution function, laser spectroscopy including LIF utilizes the (non-relativistic) Doppler shift from eigenfrequency ω 0 of target particles to absorption frequency ω: where v + v is the particle velocity parallel to the propagation direction of the incident light, and c is the speed of light. Note that the second-order Doppler effect is ignored 14 . Here, v is the contribution from the bulk fluid flow, and v is the relative velocity subject to the velocity distribution function defined in the coordinate moving at v . The frequency shift resulting from the bulk fluid flow is defined as �ω := ω 0 �v/c for convenience.
If the system is a magnetized plasma in the partial local thermodynamic equilibrium state, v of the ions is subject to the Maxwellian distribution where n i , T i , m and k B are respectively the ion density, the ion temperature, the ion mass and the Boltzmann constant. Then we obtain the spectral line shape where f (v; n i , T i )dv = f (ω; n i , T i , �ω, ω 0 )dω holds with the help of Eq. (1). In our experimental setup, we focus on the spectral line of Ar II whose eigenfrequency and mass are respectively ω 0 = 448.37 THz and m = 6.6905 × 10 −26 kg, where n i , T i , and �ω depend on space (see Methods). We assume that f (v; n i , T i ) can be observed by the LIF in any other If the fast ions at v = v b are injected into the above system, they will slow down due to the collisions of them with the thermal ions and electrons in the background plasma. Here, we suppose that the density n f of fast ions is much less than the density of thermal ions, and that v b is much greater than thermal ion velocity but much less than thermal electron velocity. Then, the velocity distribution function of fast ions will reach at the steady state, where H(· · · ) is the Heaviside step function, and v c is the crossover velocity 1,15,16 . If v < v c , the collisions with the thermal ions is dominant. If v > v c , the collisions with the thermal electrons is dominant. We also obtain the spectral line shape (1). We should mention that the line shapes of f (ω; n i , T i , �ω, ω 0 ) and g(ω; n f , ω b , ω c , �ω, ω 0 ) can be roughly the same for certain parameters and domain ( Fig. 1). While we expect that there are no fast ions in our experimental setup, it is worth challenging to identify which model represents the observed spectrum more adequately in the statistical sense, f (ω; In magnetized plasmas, the Zeeman effect is not necessarily negligible 6,8,9 . The Zeeman effect has the anisotropy dependent on the polarization of laser. When the polarization is parallel to the magnetic field, the π transition occurs. When the polarization is perpendicular to the magnetic field, the σ transition occurs. Here, we focus on the π transitions, which respectively deform f (ω; n i , T i , �ω, ω 0 ) and g(ω; n f , ω b , ω c , �ω, ω 0 ) as and where 2K is the Zeeman component number, {π k } are the transition rates, and 2δ is the Zeeman energy. Note that we ignored the dependence of �ω on ±(2k − 1)δ since (2k − 1)δ/ω 0 ≪ 1 . In our experimental setup, we www.nature.com/scientificreports/ consider that K = 3 , {π 1 , π 2 , π 3 } = {5/18, 3/18, 1/18} , and δ = 83.977 (MHz) for 0.09 T magnetic field strength 6 (see Methods).
We should emphasize that f (ω; n i , T i , �ω, ω 0 ) and f (ω; n i , T i , �ω, ω 0 , δ) are different in the internal degrees of freedom, ω or ω 0 ± (2k − 1)δ , while they reflect the same statistical behavior of ionic motion subject to If an observed spectrum is deformed, it is not always trivial to identify which factor is the origin of such a deformation, the ionic motion or the internal degrees of freedom. This ambiguity makes it more complicated to identify what factor dominates ionic motion from the observed spectrum. Under these difficulties, we statistically evaluate which form of the spectral intensity function I(ω; w) , as listed in Table , is adequate for the observed LIF spectrum, where I 0 is the background intensity. This evaluation links to identify the corresponding form of the velocity distribution function.
While many other mechanisms contribute to the spectral line broadening in LIF, such as natural broadening, power broadening, Stark broadening, and instrumental broadening, we ignore them for simplicity. This is because our experimental setup is almost the same as the setup of Ref. 6 , in which the Doppler and Zeeman effects are dominant in order. There are also many other models of ion velocity distribution function in magnetized plasma, as reviewed in Ref. 2 . Depending on the setup, they can also be candidates, not only f (v; n i , T i ) and g(v; n f , v b , v c ) , to be evaluated in terms of the statistical adequateness for the observed spectrum.
Basics of Bayesian inference. We explain the basic concept of Bayesian inference, comparing it with the conventional LIF analysis. A common approach in measuring the ion temperature and the ion flow velocity in steady-state plasmas is to find the best reproduction of a LIF spectrum by the spectral intensity function Table ) with the parameter set w = {n i , T i , �ω, I 0 } . A typical LIF spectrum of Ar II at the radial position r = 10 mm in a cylindrical helicon plasma (see Methods) and its reproduction are shown in Fig. 2a. A best-fit parameter set w for a LIF spectrum is obtained by the weighted least squares (WLS), namely minimization of the chi-squared function of w given the data set D n := {y j , ω j , σ j } n j=1 , where y j and σ j are respectively the average and the standard deviation of LIF intensity over several discharges at ω j . Now we extend this common approach to a Bayesian inference. Suppose that y 1 , y 2 , . . . , y n are independent samples taken from each conditional probability distribution In other words, if a random variable Y j is subject to p(y | ω j , σ j , w) , then it satisfies the relation Y j = I(ω j ; w) + ξ j , where ξ j is a noise subject to the Gaussian distribution N (0, σ 2 j ) . It is considered that Y 1 , Y 2 , . . . , Y n are a stochastic www.nature.com/scientificreports/ process since their realizations y 1 , y 2 , . . . , y n are originally a time-series (see Methods). Namely, ξ 1 , ξ 2 , . . . , ξ n are also a stochastic process. We expect that Eq. (9) is valid if Y j is stationary, i.e., ξ j is white, since the conditional independence of Y 1 , Y 2 , . . . , Y n , which Eq. (9) requires, is supported. In our experimental setup, we set the step size ω j+1 − ω j = 0.1 GHz for any j to sufficiently satisfy the requirement of Eq. (9) by checking the correlation length of raw data in the time domain.
In the Bayesian approach, w is also regarded as a random variable subject to the conditional probability distribution  www.nature.com/scientificreports/ where p(w) and p({y j } n j=1 | {ω j , σ j } n j=1 ) are respectively the uniform distribution as our working hypothesis and the normalizing constant. Since Eq. (10) is derived from Bayes' formula, p(w) and p(w | D n ) are respectively called the prior and posterior probability density functions.
By following the mathematical correspondence between Bayesian inference and statistical mechanics [17][18][19] , the form of p(w | D n ) is a "Boltzmann distribution" whose "inverse temperature" and "energy" for a "state" w are respectively n and χ 2 (w; D n ) as an analogy. Whereas WLS is to obtain a "ground state" as a best-fit parameter set, Bayesian inference is to obtain a whole "statistical ensemble" subject to p(w | D n ) (or to obtain p(w | D n ) itself). A "statistical ensemble", namely numerous realizations of w, subject to p(w | D n ) is shown in Fig. 2b. One can see the posterior probability density p(n i , T i | D n ) , represented by the colour gradation in the inset, tends to be higher as χ 2 (w; D n ) is smaller, where WLS solution is located around the highest density. This tendency is supported by Eq. (10); WLS solution corresponds to w that maximize p(w | D n ) . Note that this explanation is not strict but intuitive since the two values of w that respectively maximize p(w | D n ) and p(n i , T i | D n ) are almost same but strictly different 20 . In the Bayesian approach, the "ensemble" average and standard deviation are respectively taken as one of the estimators and its error bar, referred to as the posterior mean and standard deviation, since the w's "fluctuation" means the uncertainty in parameter estimation: e.g. n i = 2.0 ± 0.2 (a.u.) and T i = 0.74 ± 0.11 (eV) for Fig. 2b.

Bayesian model selection. We address the question of what form of ion velocity distribution function is
adequate to describe a observed spectrum in a statistical sense. In the Bayesian approach, a relative goodness of each I(ω; w) 's form, labeled by M l , for D n is quantified by the conditional probability where is the normalizing constant in Eq. (10), explicitly representing the dependence on M l , such as p(w | M l ) for p(w). Note that Eq. (11) is derived from Bayes' formula such that p(M) for M ∈ {M l } is the uniform distribution. We calculate p(M l | D n ) for the five different function forms listed in Table , where D n is the same spectrum shown in Fig. 2a. This result links to verify that the ion velocity distribution corresponding to the observed spectrum is a Maxwellian rather than a slowing-down. One can see from Fig. 3 that function III has the highest probability among the five forms. This probability distribution supports that the function expected to be physically valid is also statistically valid in this case. One might concern that function II (37.66%) is significantly competing with function III (38.47%), where their difference is the presence of the Zeeman effect. This competition implies that  www.nature.com/scientificreports/ the fine structure is not easy to be distinguished from the gross structure, but the observed LIF spectrum barely has the quality to do so. The posterior probability of function I significantly lower than the two above supports that there are Ar II ions at the radial position r = 10 mm. One can further confirm the uncertainty in parameter estimation. Each panel on the diagonal of Fig. 4 signifies that the marginal posterior density function of each parameter converges to Gaussian distribution. Such a convergence indicates that the amount and quality of the observed LIF spectrum are enough. Each panel on the off-diagonal of Fig. 4 signifies there is a correlation between each pair of parameters. Such a correlation comes from the form of I(ω; w) , which characterizes p(w | D n ).
Evaluation of phase-space distribution. We extend our approach to evaluate the spatial inhomogeneity in a cylindrical helicon plasma. The same analysis as above is additionally performed for each D n as observed LIF spectrum at each radial position 0 mm, 20 mm, and 30 mm (Fig. 5). One can see from Fig. 5d-f that the valid form, whose p(M l | D n ) is highest for each radial position, is function III for 0 mm, function II for 20 mm, and function I for 30 mm. The first correspond to a straightforward account of physics concerned with the experimental setup; there is the Zeeman effect induced by the external homogeneous magnetic field (Fig 3). The second indicates that the Zeeman effect is sometimes negligible, as well as other types of line broadening. One might concern a non-straightforward result that function I is valid for 30 mm. This result is a lesson that WLS solution is just an unreliable value for ion temperature and flow velocity if one straightforwardly assume function II or III for invalid situations.
We plot in Fig. 6 the radial profiles of (a) ion density n i , (b) ion temperature T i , and (c) flow velocity c�ω/ω 0 , where �ω = 0 is fixed for r = 0 mm by considering the symmetry. Two different profiles are respectively obtained under assumptions of functions II and III for comparison. They are fairly consistent, including their error bars. This consistency is connected with the fact that the posterior probabilities of functions II and III are almost the same; the Zeeman effect is sometimes negligible. One might concern significantly large errors at 30 mm, shown in Fig. 6c, corresponding to n i ≃ 0 (Fig. 6a). Such large uncertainties reflects that �ω and T i are arbitrary in the case n i ≃ 0 , namely I(ω; w) ≃ I 0 . This result is also supported by the fact that function I is valid form of the ion velocity distribution (Fig. 5f). www.nature.com/scientificreports/

Discussion
The present study sheds light on the ambiguity in identifying the velocity distribution function out of thermodynamic equilibrium states from LIF spectra. We demonstrate that our approach can resolve the ambiguity. The success of our approach relies on reliable LIF measurements. Our approach is also applicable to other spectroscopic techniques, such as absorption spectroscopy and two photons absorption LIF 14 .
Our approach may not work as expected if the signal-to-noise ratio (SNR) is not sufficiently high. The uncertainty in parameters and function forms is connected with the SNR in the Bayesian inference. It is pointed out that there are thresholds of the SNR permissible for parameter estimation 21 and model selection 22 in similar contexts. The systematic error from experimental artefacts, e.g. the thermal drift of laser and the reproducibility of plasma production, may also bring unexpected results since we suppose that there are only random errors in our formulation (Eq. (9)). The plasma instabilities and fluctuation may be sources of error depending on the setup. In order to obtain more reliable results, it is essential to minimize these systematic errors and disturbances.

PANTA experiments.
The experiment was carried out in a linear magnetized plasma device, PANTA 23 .
The vacuum vessel has the length of 4050 mm and the diameter of 450 mm. 0.09 T of homogeneous magnetic field is controlled by 17 pairs of Helmholtz coils. Argon plasma with 50 mm radius is produced by helicon wave with 3 kW at 7 MHz radio frequency, and neutral gas pressure is fixed at 0.1 Pa, where the streamer structure are www.nature.com/scientificreports/ observed 24 . The discharge duration is 2.0 s. Typical central electron density and temperature are ∼ 1 × 10 19 m −3 and 3 eV, respectively 25 .

LIF measurements.
Here, we briefly review the LIF measurement system in PANTA. A tunable diode laser tuned at around 668.61 nm is injected perpendicular to the magnetic fields and utilized for exciting the 3d 4 F 7/2 level to the 4p 4 D 5/2 level in Ar II in the experiments 5 . The laser wavenumber is swept around 668.61 ± 0.01 nm by triangle waveform with 0.1 Hz of modulation frequency. The effect of laser thermal drift has been calibrated using a Fabry-Pérot interferometer and an iodine cell. The LIF signal is collected by a photo multiplier tube and amplify and averaging by lock-in amplifier. Since the error of the LIF signal is large, we performed 180 discharges at each radial position in the experiments The LIF signals at 0.2 s after the start of discharge, when the turbulence becomes quasi-steady state, is used to evaluate the ion velocity distribution function. For more details of the LIF system, see Ref. 9 .  26,27 , where the total MC sweeps were 10,000 after the burn-in. To calculate p({y j } n j=1 | {ω j , σ j } n j=1 ) , we utilized the bridge sampling 28,29 .