Collective optical Thomson scattering (TS) is a diagnostic commonly used to characterize plasma parameters. These parameters are typically extracted by a fitting algorithm that minimizes the difference between a measured scattered spectrum and an analytic spectrum calculated from the velocity distribution function (VDF) of the plasma. However, most existing TS analysis algorithms assume that the VDFs are Maxwellian, and applying an algorithm that makes this assumption does not accurately extract the plasma parameters of a non-Maxwellian plasma due to the effect of non-Maxwellian deviations on the TS spectra. We present new open-source numerical tools for forward modeling analytic spectra from arbitrary VDFs and show that these tools are able to more accurately extract plasma parameters from synthetic TS spectra generated by non-Maxwellian VDFs compared to standard TS algorithms. Estimated posterior probability distributions of fits to synthetic spectra for a variety of example non-Maxwellian VDFs are used to determine uncertainties in the extracted plasma parameters and show that correlations between parameters can significantly affect the accuracy of fits in plasmas with non-Maxwellian VDFs.

## I. INTRODUCTION

Thomson scattering (TS) refers to the scattering of electromagnetic radiation by a collection of many charged particles, such as a plasma.^{1,2} Optical TS, in which the probing radiation consists of optical wavelengths, is a widely used, *in situ*, non-perturbative diagnostic tool for characterizing plasmas, with applications ranging from laboratory astrophysics^{3–7} to fusion plasmas.^{8–12}

In collective TS, the parameters associated with the plasmas of interest cause incident optical radiation to be scattered by electron plasma waves (EPWs) and ion acoustic waves (IAWs). The fluctuations in the charge density produce TS spectral features,^{1} which can be measured and analyzed in order to extract information about the particle velocity distribution functions (VDFs). Most TS analysis tools assume that the plasma is thermalized so that the electron and ion VDFs are both Maxwellian with the corresponding parameters, such as the temperatures, *T*_{e} and *T*_{j}, and the densities, *n*_{e} and *n*_{j}, of the electron and ion populations. These Maxwellian parameters can then be inferred from the spectrum by an algorithm that minimizes the difference between the measured scattered spectrum and an analytic spectrum generated from these Maxwellian VDFs.

Non-Maxwellian VDFs have directly been observed with TS diagnostics in recent experiments on high-energy-density (HED) plasmas^{13–16} and laser-driven collisionless shocks.^{4,17} In addition, non-Maxwellian VDFs are expected to occur in the divertor region of tokamak devices.^{18} The deviations of the VDF from a Maxwellian can be crucial to understanding the plasma dynamics. Previous theoretical studies have discussed the form of the TS spectra generated from non-Maxwellian distributions, such as a two-stream distribution^{19,20} and a super-Gaussian distribution,^{21} as well as how the resulting spectra deviate from their Maxwellian counterparts. Milder *et al.*^{22} also showed that fitting the TS spectrum produced by a super-Gaussian electron VDF with a Maxwellian model can give inaccurate fitted plasma parameters. The inaccurate fitting is due to how the non-Maxwellian deviations affect the strength of Landau damping at the location of the Thomson spectral peaks. Consequently, existing tools designed for Maxwellian VDFs are insufficient to analyze the TS spectra from plasmas with non-Maxwellian VDFs.

In this paper, we investigate how non-Maxwellian VDFs impact TS spectra with the aid of two new numerical tools that we have developed.^{23} The first tool is a “forward model,” which computes the analytic TS spectrum from a set of arbitrary discretized particle VDFs, and the second tool is a “fitting algorithm,” which extracts non-Maxwellian plasma parameters from a TS spectrum by iteratively applying the non-Maxwellian forward model at different points in parameter space. These open-source tools expand on the work of Milder *et al.* by enabling the study of arbitrary non-Maxwellian VDFs.

For the fitting algorithm, we examine the use of two schemes to optimize the parameter space exploration. One scheme is differential evolution (DE), which attempts to continuously optimize a candidate solution by combining previous solutions. This scheme can search a large parameter space and is often computationally more efficient than “brute force” methods, but it is susceptible to uncertainties caused by correlations between parameters in the solution. The second scheme is a Monte Carlo Markov Chain (MCMC). While generally computationally more expensive than DE, MCMC is well-suited for exploring very large parameter spaces that can be used to estimate the uncertainties in the best-fit parameters from the DE scheme.

In Sec. II, we review the TS theory that forms the basis for our method. The numerical tools we developed to implement this routine are discussed in Secs. III and IV. We also describe a process for testing our method and comparing it with an open-source method that assumes Maxwellian VDFs, of which the results are presented in Sec. V. Finally, in Sec. V B, we analyze the uncertainty and confidence associated with our fitting algorithm and discuss possible areas of improvement. Our conclusions are summarized in Sec. VI.

## II. THOMSON SCATTERING THEORY

*et al.*

^{1}TS occurs when the absorption of an incident photon causes a charged particle to undergo acceleration, which then induces Larmor radiation as the charge emits a photon. In the charge’s rest frame (the primed frame), the frequency $\omega s\u2032$ of the scattered photon is equal to the frequency $\omega i\u2032$ of the incident photon. In the (unprimed) lab frame,

*ω*

_{s}can be solved for by computing the primed frame solution and applying the appropriate Doppler shifts, which are a function of the particle velocity

**v**as seen in the lab frame, as well as the incident and scattered wavevectors

**k**

_{i}and

**k**

_{s}, respectively. If we define

*ω*≡

*ω*

_{s}−

*ω*

_{i}and

**k**≡

**k**

_{s}−

**k**

_{i}, the final result is shown to be

*P*has the following proportionality:

*S*(

**k**,

*ω*) is the spectral density function, which contains the dependence on the velocity distributions of the charged particles in the scattering volume. Note that when discussing TS forward models and fitting algorithms, the “TS spectrum” being computed and fitted is either the spectral density function itself or the scattered power as in Eq. (2), depending on the data being fit. For the purposes of this paper, “TS spectrum” refers to the normalized scattered power unless otherwise specified. In general, the spectral density function can be written in terms of the normalized VDFs of the electrons and ions in the plasma,

*f*

_{eo}is the one-dimensional projection of the electron VDF in the direction of measurement,

*f*

_{jo}are the one-dimensional projections of the ion VDFs, with

*j*indexing the ion species,

*Z*

_{j}and

*n*

_{j0}are the charge and density of ion species, respectively, and

*N*is the combined density of all ions. The electron susceptibility

*χ*

_{e}and the ion susceptibilities

*χ*

_{j}are functions of

**k**and

*ω*given by

*n*

_{e0}is the electron density,

*e*is the electric charge,

*m*

_{e}is the electron mass, and the integrals can be performed along a Landau contour, which deviates from the real axis just enough to avoid the pole at

*v*=

*ω*/

*k*. The longitudinal dielectric function

*ϵ*is given by

*ω*

_{s}≈

*ω*

_{i}. Therefore,

**k**

_{s}and also

**k**will not vary significantly. If we make the approximation that the direction of

**k**is effectively fixed, then

*χ*

_{e}and

*χ*

_{j}will only be sensitive to the 1D projected VDF in the direction of

**k**and Eqs. (4) and (5) can be rewritten as integrals over scalars.

*α*, defined in terms of the electron Debye length

*λ*

_{De}as

*α*≪ 1, then

*λ*≪

*λ*

_{De}and the incident radiation effectively sees the electrons as free (unbound), leading to the non-collective regime. If

*α*≳ 1, then the radiation sees the Debye-shielded charges and, therefore, the correlations between the motions of the electrons, leading to collective scattering. The collective TS spectrum can be broken down into two regimes. At high |

*ω*|, the heavier ions are unable to respond, so the electron component as labeled in Eq. (3) dominates. At low |

*ω*|, both the electrons and ions are able to respond, but the ion component tends to dominate. The electron component is generally still non-trivial, which becomes important when applying a fitting algorithm. We will refer to the high |

*ω*| regime as the EPW spectrum, as it is dominated by scattering from EPWs, and the low |

*ω*| regime as the IAW spectrum, as it is dominated by scattering from IAWs. Due to the difference in frequency scales between the electron- and ion-dominated components, the EPW and IAW spectra are usually measured with separate spectrometers in experiments,

^{4,24}so treating them separately is justified. In addition, the measurement of the EPW spectrum usually involves the placement of a notch filter to block out the low-

*ω*portion of the spectrum as it is much brighter than the high-

*ω*electron contributions. Our code takes this notch filter into account.

Because Eq. (3) can be used to compute the TS spectrum from an arbitrary set of VDFs, in principle, it may be possible to invert this process and infer the VDFs from a TS spectrum, assuming that the TS spectrum is a non-degenerate function of the VDFs. In practice, this inversion can be unreliable even if the TS spectrum is nearly degenerate or if the VDFs are described by too many parameters.

This inversion process is known to work if the VDFs in question are Maxwellian,^{1} and most tools developed for the TS analysis use this Maxwellian assumption. However, these tools cannot accurately reconstruct the VDFs when this assumption is violated, as described in the following sections.

### A. Example VDFs

Although in theory, it is possible to define a TS spectrum using arbitrary VDFs, the examples presented in the figures in this paper focus on a representative selection of non-Maxwellian VDF models relevant to plasma physics. These VDF models as well as their plasma parameters are discussed below. Because the TS spectrum is only sensitive to the one-dimensional projections of the VDFs in the direction of **k**, the one-dimensional projections corresponding to the example VDF models are also presented below.

*k*

_{B}is the Boltzmann constant,

*T*is the temperature, and

**v**

_{D}is the bulk drift velocity.

**k**to lie along the

*x*axis, the relevant one-dimensional projection is

*v*

_{Dx}is the bulk drift velocity in the

**k**direction.

A VDF can also be composed of a linear combination of two or more Maxwellians, which have different temperatures and/or drift velocities.

*m*. Here,

*v*

_{D}is the drift velocity, Γ() is the Gamma function, and the spectral index

*κ*is a measure of the non-Maxwellian deviation.

^{25}Note that the usual notion of temperature does not apply to non-Maxwellian distributions as they are not thermalized, but we can still define an equivalent temperature

*T*in terms of the variance of the VDF,

*v*

_{D}is defined for arbitrary VDFs as

*κ*→

*∞*, the kappa distribution approaches a Maxwellian with the same drift velocity and equivalent temperature, while at

*κ*= 3/2, the distribution collapses and becomes undefined. Taking a linear combination of a hot kappa distribution and a cooler Maxwellian results in a core-halo distribution, which is commonly observed in the solar wind.

^{26}

^{27}

^{22}The general form looks similar to a Maxwellian or Gaussian distribution but with the power in the exponent as an additional free parameter. The super-Gaussian distribution takes the form

^{22}

**v**

_{D}and additional parameters

*v*

_{p}and

*p*. The characteristic velocity

*v*

_{p}is related to the equivalent temperature by

*p*= 2, this reduces to a Maxwellian [Eq. (8)] and

*v*

_{p}becomes the thermal speed $vp=vth\u22612kBT/m$. For

*p*≠ 2, the parameter

*v*

_{p}remains linearly related to the thermal speed but multiplied by a constant factor, which depends on

*p*. The one-dimensional projection of the super-Gaussian can be written in an integral form,

### B. *χ* for different VDF models

*χ*functions defined in Eqs. (4) and (5), note that if we perform a change of variables $u=v/vth2$ and $\xi =\omega /kvth2$ and operate in the regime where

*k*is nearly fixed, then we can write

*χ*as a function of

*ξ*, up to some constant factors determined by the plasma parameters,

These dimensionless integrals have the temperature dependence factored out, so they only depend on the overall form of the VDF and will differ for different VDF models. Figure 1 shows the dimensionless *χ*(*ξ*) integrals for three different VDF models for reference. It is helpful to note that the real part of *χ* represents the dispersion of an electrostatic wave at given *ω* through the plasma, while the imaginary part represents the Landau damping on electrostatic waves.

The *χ* functions significantly impact the features of the TS spectrum. For instance, a larger imaginary part of *χ*, corresponding to stronger Landau damping, is associated with broader spectral peaks.^{22} This can be seen in Fig. 2, which shows several analytic TS spectra from Maxwellian and drifting Maxwellian VDFs. The location of the EPW spectral peaks depends on the dispersion relation of EPWs,^{1} which only weakly depends on the temperature at optical wavelengths (assuming that the charge density is fixed). Therefore, increasing the temperature leaves the location of the spectral peaks essentially unchanged in velocity space but increases the thermal velocity and, therefore, moves the spectral peaks to lower *ξ*, which corresponds to higher imaginary *χ* (for *ξ* ≳ 1). Thus, we expect broader spectral peaks at higher temperatures, which are seen in the EPW TS spectra in Fig. 2.

Non-Maxwellian deviations in the VDFs impact *χ* even if the temperature and density are held constant, which then affect the TS spectrum. A purely Maxwellian fitting algorithm can only vary the Maxwellian temperature and density, so when fitting a TS spectrum with non-Maxwellian deviations in *χ*, the algorithm will attempt to compensate for those non-Maxwellian deviations by changing the temperature and density in order to affect *χ*. This leads to incorrect fitting of the plasma parameters. This can be illustrated following the above example with temperature. For a given *ξ* ≳ 1, imaginary *χ* (Landau damping) is smaller for a super-Gaussian compared to a Maxwellian distribution, as shown in Fig. 1. Consequently, the same temperature (as defined in Sec. II A) will lead to narrower spectral peaks for the super-Gaussian.

## III. ARBITRARY FORWARD MODEL

Here, we describe a forward model that maps a given set of discretized arbitrary VDFs to a corresponding TS spectrum, assuming that the VDFs are in quasi-equilibrium and that the plasma is unmagnetized.^{1} The plasma physics scientific Python package PlasmaPy^{28} includes a built-in function that computes the TS spectrum given a set of Maxwellian plasma parameters (the “Maxwellian forward model”), which includes ion and electron temperatures, densities, and drift velocities. We expand upon this code and develop a numerical function that computes the TS spectra for arbitrary VDFs, which we refer to as the “arbitrary forward model.” For a given VDF, the arbitrary forward model accepts an array of velocity values {*v*_{j}} and the corresponding VDF values {*f*_{j}}, which are meant to represent *f*_{j} = *f*(*v*_{j}), where *f*(*v*) is the true continuous VDF. The Maxwellian forward model in PlasmaPy was used as a benchmark for testing and confirming the validity of our arbitrary forward model.

### A. Numerical implementation

*v*=

*ω*/

*k*. We note that the susceptibility integrals can be recast into the form

*v*→

*u*(

*v*). This integral can be rewritten using Plemelj’s formula

^{1}as

*∫*refers to the Cauchy principal value of the integral, defined as

*ϵ*=

*ϕ*, we can calculate the principal value by standard numerical integration in the ranges [

*a*,

*ξ*−

*ϕ*] and [

*ξ*+

*ϕ*,

*b*], plus a correction term,

*g*(

*ξ*) is proportional to

*∂f*/

*∂v*, which is numerically computed from the discretized input VDFs using a finite difference scheme to fourth order precision. The two definite integrals in Eq. (24) cross no poles, so they can be evaluated using a standard numerical integration scheme, such as a Riemann sum. To minimize the number of integration points while maintaining precision, we implement a scheme in which the points are finely spaced close to the pole, where the integrand varies rapidly, but coarsely spaced away from the pole, where the integrand varies slowly.

### B. Numerical errors

The estimation of the *χ* integrals as performed in Eq. (24) as part of the arbitrary forward model introduces some numerical errors into the TS spectrum. Our algorithm uses the range of velocities defined for the input VDFs to determine an integration range. It also obtains the values of the integration points by interpolating from the input VDFs. This means that appropriate resolution of the input VDFs is necessary to minimize the error of the forward model. This constrains both the array spacing *δv* of the input VDFs and the total range over which the input VDFs are defined.

For a VDF *f*(*v*), the array spacing must be sufficiently small to resolve the features of the VDF and to precisely compute *∂f*/*∂v*. In general, *δv* should be less than a characteristic velocity associated with the resolution of *f*(*v*), found by taking the ratio of the VDF to its derivative: *δv* ≤ *f*(*v*)/*f*′(*v*). The total range of the velocity array must be large enough so that the forward model includes all important features of the VDF. This can be expressed as *f*(*v*)/max(*f*) ≪ 1 for all *v* outside the array range, where max(*f*) is the maximum value that *f*(*v*) takes over all *v*. Both of these conditions are necessary for the numerical calculation of the spectral density. If the input VDF spacing is too sparse, then the *χ* integration may not have enough points to be accurate. If the VDF range is too small, then the finite integration is not a good approximation of the integral from −*∞* to *∞*. Better quantification of the effects of the input velocity array on numerical integration will be the subject of future studies.

We quantify the error in the Maxwellian case by comparing the output spectra generated by the arbitrary forward model using Maxwellian input VDFs with the corresponding output spectra from the Maxwellian forward model. Figure 3 shows the comparison of the forward-modeled spectra for examples with VDFs that are Maxwellians or linear combinations of Maxwellians. There are a few small differences in the spectra, but they have good agreement overall, validating the arbitrary forward model. The pointwise relative error of the forward-modeled spectra from the arbitrary forward model was ≲5% for the examples shown in Fig. 3. In addition, those differences can be decreased arbitrarily by increasing the resolution of the *χ* integration scheme.

To illustrate the numerical effects of velocity resolution and VDF range, Fig. 4 shows the output spectra when the arbitrary forward model is applied to the same Maxwellian VDFs (example 1 in Table I), but where the velocity arrays on which the VDFs are defined have been modified. The ranges of the velocity arrays have been changed while leaving the array length fixed in both cases. In example 5, the velocity range is much larger and the VDF is effectively defined by very few velocity bins. We refer to this as the “narrow” case. In example 6, the velocity range is smaller so that the tails of the VDFs are cut off. We refer to this as the “wide” case. In the narrow case, the array spacing *δv* is several orders of magnitude larger than the characteristic velocity *f*(*v*)/*f*′(*v*) at some values of *v*, so the VDF is not well-resolved. In the wide case, the VDFs are wider than the velocity array. The VDFs are cut off at a factor of $\u223c0.01$ of their maxima, which is much larger compared to ≲10^{−10} in the well-resolved case shown in Fig. 3. In both cases, the resulting computed spectra deviate significantly from the Maxwellian forward model.

Ex. # . | T_{e} (eV)
. | v_{e} (m/s)
. | T_{i} (eV)
. | v_{i} (m/s)
. | n_{e} (cm^{−3})
. | n_{i} (cm^{−3})
. | Resolution (δv_{e}, δv_{i}) (m/s)
. |
---|---|---|---|---|---|---|---|

1 | 300 | 0 | 100 | 0 | 4 × 10^{18} | 4 × 10^{18} | (2 × 10^{5}, 8 × 10^{3}) |

2 | 100 | −3 × 10^{6} | 150 | 0 | 4 × 10^{18} | 4 × 10^{18} | ⋯ |

3 | 200 | 0 | 50 | 2 × 10^{5} | 4 × 10^{18} | 4 × 10^{18} | ⋯ |

4 | 300 | 0 | (100, 200) | (0, 3) × 10^{5} | 4 × 10^{18} | (1.6, 2.4) × 10^{18} | (2 × 10^{5}, 8 × 10^{3}) |

5 | 300 | 0 | 100 | 0 | 4 × 10^{18} | 4 × 10^{18} | (8 × 10^{6}, 4 × 10^{5}) |

6 | 300 | 0 | 100 | 0 | 4 × 10^{18} | 4 × 10^{18} | (8 × 10^{4}, 1.2 × 10^{3}) |

Ex. # . | T_{e} (eV)
. | v_{e} (m/s)
. | T_{i} (eV)
. | v_{i} (m/s)
. | n_{e} (cm^{−3})
. | n_{i} (cm^{−3})
. | Resolution (δv_{e}, δv_{i}) (m/s)
. |
---|---|---|---|---|---|---|---|

1 | 300 | 0 | 100 | 0 | 4 × 10^{18} | 4 × 10^{18} | (2 × 10^{5}, 8 × 10^{3}) |

2 | 100 | −3 × 10^{6} | 150 | 0 | 4 × 10^{18} | 4 × 10^{18} | ⋯ |

3 | 200 | 0 | 50 | 2 × 10^{5} | 4 × 10^{18} | 4 × 10^{18} | ⋯ |

4 | 300 | 0 | (100, 200) | (0, 3) × 10^{5} | 4 × 10^{18} | (1.6, 2.4) × 10^{18} | (2 × 10^{5}, 8 × 10^{3}) |

5 | 300 | 0 | 100 | 0 | 4 × 10^{18} | 4 × 10^{18} | (8 × 10^{6}, 4 × 10^{5}) |

6 | 300 | 0 | 100 | 0 | 4 × 10^{18} | 4 × 10^{18} | (8 × 10^{4}, 1.2 × 10^{3}) |

The pointwise relative error of the well-resolved, narrow, and wide cases is shown in Fig. 5. Due to the *χ* integration over the entire VDF, every point on the TS spectrum is affected by the entire VDF. Even in the wide case, where a subset of the VDF is very well-resolved, the pointwise error of the arbitrary forward model is worse across the entire TS spectrum.

## IV. FITTING ALGORITHMS

The current release of the plasma physics Python package PlasmaPy^{28} contains a built-in function, which can fit the Maxwellian plasma parameters to the TS spectra using the DE algorithm from the Python package lmfit.^{29} This algorithm can also fit to sums of Maxwellians. Because the computation of the susceptibility function for Maxwellians is optimized by tabulation,^{30} using the Maxwellian forward model makes the fitting significantly faster than it would be with the arbitrary forward model. Therefore, the main use case of the arbitrary forward model for fitting would be in fitting of TS spectra from non-Maxwellian VDFs.

There are several possible approaches to a non-Maxwellian iterative fitting algorithm, which produce fits to varying degrees of arbitrariness. If we had sufficient computing power, we could, in theory, treat each point on a discretized VDF array as its own parameter and then fit it to that set of parameters. However, implementing this approach with a sufficient resolution of the velocity axis becomes computationally infeasible due to the large number of free parameters. Another potential approach would be to provide a number of pre-defined parameterized non-Maxwellian VDF models with low-dimension parameter spaces that the algorithm can fit to. This has the issue of being too restrictive, and a better method in this case would be to design individual forward models based around each of the VDF models.

The approach our algorithm employs is to accept custom user-defined parameterized VDF models and then to fit the input TS spectra assuming that the VDFs are in the form of those user-defined models. This allows the user to fit to different non-Maxwellian models while still keeping the parameter space small. We compared the fits obtained by this approach with the fits obtained from assuming Maxwellian VDFs to determine how well the arbitrary method performs.

As discussed in Sec. II, the EPW spectrum is negligibly affected by the ion VDFs, so the EPW spectrum is first used to fit the electron parameters while holding the ion parameters fixed at some arbitrary values that would be reasonable for the system. The IAW spectrum depends on both electron and ion parameters, so to reduce the time needed for fitting, the electron parameters can be fixed at the values obtained from fitting the EPW spectrum and the remaining ion parameters are fitted using the IAW spectrum. If necessary, this process could be further iterated with additional parameter constraints to refine the fit as needed.

### A. Procedure for comparing fitting algorithms

We apply a general procedure to test both fitting algorithms and compare the results. First, we prepare synthetic electron and ion VDFs, which have the same form as physically relevant VDFs. We then compute the EPW and IAW TS spectra of a plasma with these VDFs. Gaussian noise is added to make the resulting spectra more realistic, with a standard deviation given by *σ*_{n} = 0.1 max(*P*(*λ*)), with *P*(*λ*) being the TS spectrum. The VDFs are then reconstructed based on the fitted parameters and are compared with the initial input VDFs, while the parameters and their variances are compared with the input parameters. After obtaining the results from the fitting algorithms, an MCMC sampler is used to explore the parameter space and compute the one- and two-dimensional posterior probability distributions (PPDs) of the parameters in order to estimate the confidence in the fitted values.

After fitting the synthetic spectra, we analyze and compare the accuracy of the fits by calculating the *χ*^{2} statistic between the best-fit VDFs and the input VDFs, as well as the percent error of quantities, such as density, bulk flow velocity, and the equivalent temperature of the VDFs. We characterize the error bars associated with the fit and use these to search for correlations and degeneracies in the arbitrary forward model using an MCMC sampler.^{31}

### B. Validating the arbitrary fitting algorithm

*χ*

^{2}goodness of fit,

*I*and fitted function

*F*. In the case of the TS spectra, the Gaussian noise that is artificially added to the input spectra will increase

*χ*

^{2}on its own, so the usual definition of

*χ*

^{2}is normalized by $\sigma n2$ to account for this and allow for comparison between the different fitting examples. Therefore, for the TS spectra, we use

*χ*

^{2}< 1, if the randomly generated synthetic noise happens to contribute less than the expected squared-error of $\sigma n2$. In Table III, we see that the Maxwellian and arbitrary fitting algorithms are approximately on a par with each other when fitting Maxwellians.

Ex. # . | Ion species . | T_{e} (eV)
. | v_{e} (m/s)
. | T_{i} (eV)
. | v_{i} (m/s)
. | n_{e} (cm^{−3})
. | n_{i} (cm^{−3})
. | κ_{e}
. | p_{C}
. |
---|---|---|---|---|---|---|---|---|---|

7 | p | 200 | 10^{6} | 50 | 10^{5} | 4 × 10^{18} | 4 × 10^{18} | ⋯ | ⋯ |

8 | p | 200 | 0 | 50 | 10^{5} | 4 × 10^{18} | 4 × 10^{18} | 2 | ⋯ |

9 | (p, $C125+$) | 150 | 0 | (100, 200) | (2 × 10^{5}, 0) | 4 × 10^{18} | (2.4, 0.267) × 10^{18} | ⋯ | 3.5 |

Ex. # . | Ion species . | T_{e} (eV)
. | v_{e} (m/s)
. | T_{i} (eV)
. | v_{i} (m/s)
. | n_{e} (cm^{−3})
. | n_{i} (cm^{−3})
. | κ_{e}
. | p_{C}
. |
---|---|---|---|---|---|---|---|---|---|

7 | p | 200 | 10^{6} | 50 | 10^{5} | 4 × 10^{18} | 4 × 10^{18} | ⋯ | ⋯ |

8 | p | 200 | 0 | 50 | 10^{5} | 4 × 10^{18} | 4 × 10^{18} | 2 | ⋯ |

9 | (p, $C125+$) | 150 | 0 | (100, 200) | (2 × 10^{5}, 0) | 4 × 10^{18} | (2.4, 0.267) × 10^{18} | ⋯ | 3.5 |

. | Maxwellian fitting algorithm . | Arbitrary fitting algorithm . | ||||||
---|---|---|---|---|---|---|---|---|

Ex. # . | $\chi EPW2$ . | $\chi eVDF2$ . | $\chi IAW2$ . | $\chi iVDF2$ . | $\chi EPW2$ . | $\chi eVDF2$ . | $\chi IAW2$ . | $\chi iVDF2$ . |

7 | 1.005 | 2.97 × 10^{−20} | 1.203 | 4.87 × 10^{−16} | 1.006 | 3.58 × 10^{−20} | 1.200 | 6.19 × 10^{−16} |

8 | 1.534 | 2.42 × 10^{−17} | 1.507 | 3.265 × 10^{−15} | 1.092 | 8.38 × 10^{−19} | 1.022 | 1.403 × 10^{−16} |

9 | 1.147 | × 10^{−19} | 1.010 | (0.171, 1.72) × 10^{−12} | 1.148 | 9.01 × 10^{−20} | 0.956 | (2.399, 0.678) × 10^{−14} |

. | Maxwellian fitting algorithm . | Arbitrary fitting algorithm . | ||||||
---|---|---|---|---|---|---|---|---|

Ex. # . | $\chi EPW2$ . | $\chi eVDF2$ . | $\chi IAW2$ . | $\chi iVDF2$ . | $\chi EPW2$ . | $\chi eVDF2$ . | $\chi IAW2$ . | $\chi iVDF2$ . |

7 | 1.005 | 2.97 × 10^{−20} | 1.203 | 4.87 × 10^{−16} | 1.006 | 3.58 × 10^{−20} | 1.200 | 6.19 × 10^{−16} |

8 | 1.534 | 2.42 × 10^{−17} | 1.507 | 3.265 × 10^{−15} | 1.092 | 8.38 × 10^{−19} | 1.022 | 1.403 × 10^{−16} |

9 | 1.147 | × 10^{−19} | 1.010 | (0.171, 1.72) × 10^{−12} | 1.148 | 9.01 × 10^{−20} | 0.956 | (2.399, 0.678) × 10^{−14} |

The fitted electron density, electron temperature, and ion temperature are all within 5% of the true values for both algorithms, further showing that both algorithms can accurately fit the TS spectra from Maxwellians. The fitted electron parameters from both algorithms and their values relative to the true values are also graphically shown in Fig. 9.

## V. RESULTS

### A. Comparisons between best-fit VDFs

After validating the arbitrary fitting algorithm, we test it on synthetic TS spectra from plasmas with non-Maxwellian VDFs, which are computed using the arbitrary forward model. The resulting fits are compared with those generated by applying the Maxwellian fitting algorithm to the same data. The plasma parameters for these non-Maxwellian parameters are also given in Table II.

First, we study the combination of a non-thermal kappa distribution for the electron VDF and a drifting Maxwellian distribution for the ion VDF. As can be seen in Fig. 7, the electron VDF is not well-fit by a Maxwellian.

While the ions are Maxwellian, the ion parameters are still not correctly fit to the Maxwellian fitting algorithm. This is because the IAW spectrum has a non-trivial dependence on the electron VDF so that the errors in the electron VDF effectively propagate into the IAW fitting. We also see this in Table III, as the $\chi IAW2$ and $\chi iV DF2$ values are significantly higher for the Maxwellian fitting algorithm than for the arbitrary fitting algorithm. For the arbitrary fitting algorithm, because the correct models are used for both the electrons and the ions, the VDFs are both well-fit. However, if the wrong VDF model were to be used for the electrons, the arbitrary fitting algorithm would suffer from the same error propagation issue. This emphasizes the importance of using the correct VDF models for fitting.

In addition, the plasma parameters are not accurately recovered by the Maxwellian fitting algorithm. For instance, the fitted electron density has an error of ~11% and the fitted electron temperature has an error of ~35%, which is much more than the error in the fitted plasma parameters when fitting TS from Maxwellians. In comparison, the fitted electron temperature and density from the arbitrary fitting algorithm are within 5% error, which is a significant improvement over the Maxwellian algorithm. The fitted electron parameters from both algorithms and their values relative to the true values for this example are shown in Fig. 10.

Second, we attempt to fit the synthetic spectrum of a drifting Maxwellian electron VDF, a drifting Maxwellian proton VDF, and a super-Gaussian carbon ion (C^{+5}) VDF. The relative densities of the protons and carbon ions are chosen to make the plasma quasi-neutral. In Fig. 8, we see that the electron VDF is well-fit by the Maxwellian fitting algorithm, but the proton and carbon VDFs are fitted inaccurately. Similar to the previous example where the poorly fit non-Maxwellian electron VDF impacted the proton VDF fitting, here, the non-Maxwellian carbon VDF causes a similar effect.

Again, there is a discrepancy between the fitting algorithms with the relative errors of the fitted plasma parameters. Here, the Maxwellian fitting algorithm underestimates the proton temperature by over 50%, while the arbitrary fitting algorithm overestimates it by ~23%. The larger errors are to be expected due to the greater number of parameters being fitted, but we see that the arbitrary fitting algorithm is still significantly more accurate than the Maxwellian algorithm.

Although the best-fit ion VDFs for the Maxwellian fitting algorithm are quite different from the corresponding best-fit VDFs from the arbitrary fitting algorithm, we see that the resulting best-fit IAW spectra look sufficiently similar to be within noise effects of each other. This can also be seen in Table III, as the $\chi IAW2$ values are similar for the two fitting algorithms, but the $\chi iV DF2$ values are significantly higher for the Maxwellian fitting algorithm. These results suggest a near-degeneracy in the TS forward model.

The existence of near-degeneracies in the TS forward model raises the risk of getting a good fit on the TS spectrum, which corresponds to an inaccurate VDF. This can be mitigated by restricting the fitting algorithm to a physically motivated VDF model and using other diagnostic results to put limits on the range of plasma parameters. Section V B discusses how to identify degeneracies arising from different VDF models.

### B. Uncertainty analysis

After finding the best-fit plasma parameters, we analyze the robustness and uncertainty in the fits. This allows us to determine which aspects of the fits we can be confident in and to reveal TS spectrum degeneracies that could make the fits inaccurate. The Python module emcee^{31} was used to estimate the posterior probability distributions (PPDs) of the fits using an MCMC sampling algorithm, which is given initial values at the best-fit values from the arbitrary fitting algorithm. Using the PPDs, we estimate error bars for each fitted parameter and show that when fitting TS spectra from non-Maxwellian distributions, fits using the arbitrary fitting algorithm derived plasma parameters with greater accuracy than the Maxwellian fitting algorithm.

Figure 9 shows the PPD associated with the fitting of the EPW spectrum in Fig. 6. We see that the 1D PPD projections [Figs. 9(a), 9(c), and 9(f)] approximately follow Gaussian distributions. The two-dimensional slices [Figs. 9(b), 9(d), and 9(e)] also look approximately symmetric, indicating that the variables are effectively fitted independently of each other. From this, we conclude that the fitted parameters that were output by the arbitrary fitting algorithm are accurate, as there is one clear location in parameter space where the posterior probability density reaches a maximum, and all marginal distributions also reach their maxima. The 95% confidence intervals (dashed lines) for each parameter are shown for reference.

Both fitting algorithms produce approximately equal-sized uncertainties for each parameter, and the true values are contained within these uncertainties. This benchmarks the arbitrary fitting algorithm to be approximately on a par with the Maxwellian fitting algorithm when fitting TS spectra from Maxwellians. Some of the best-fit values in Fig. 9 visually appear far away from the true parameters, but this is due to a rescaling of the relevant axes in the figure in order to see the PPD features and does not represent the actual ranges of parameters over which the algorithm was allowed to explore (this is especially the case with the velocity parameter, which was fit over a large range). For instance, the fitting algorithm used to produce the results of Fig. 9 ranged over electron temperatures of 10–2000 eV, drift velocities of −10^{7} to 10^{7} m/s, and electron densities of 10^{17}–10^{19} cm^{−3}. When taking the full ranges into account, the errors of the best fit parameters are on the order of $\u223c1%$ of the total parameter range, which shows that the algorithms find the correct plasma parameters even when allowed to explore over a very large region of parameter space.

Figure 10 shows the PPD corresponding to the fitting of the electron parameters from the EPW spectrum in Fig. 7 using the arbitrary fitting algorithm. Unlike the Maxwellian case, the 1D marginal distributions [Figs. 10(a), 10(c), and 10(f)] are no longer Gaussian and the two-dimensional slices [Figs. 10(b), 10(d), and 10(e)] illustrate a strong correlation between a number of fitted plasma parameters, especially the electron temperature and the spectral index *κ*. This is because both the temperature and *κ* affect the width of the VDF, so they are degenerate. We see that the 95% confidence intervals for *T* and *κ* are large, which, in the absence of other information, might suggest that the fit is not robust. However, from the full PPD, we see that the single-parameter confidence intervals are misleading due to the correlated parameters, as the PPD remains concentrated in a relatively small region of the full parameter space. It is possible that if the initial guesses for the DE fitting algorithm are chosen poorly, the algorithm would land somewhere along the curve where the value of the PPD is high but far from the true values, thereby resulting in inaccurate fitted parameters. In practice, physical constraints obtained by analyzing other experimental data must be applied as boundaries on the fitting parameters to break this degeneracy.

We also see that the Maxwellian-fitted electron temperature is significantly lower than the true equivalent temperature of the kappa distribution. This is due to the fact that at the location of the EPW spectral peaks, corresponding to *ξ* ≈ 1.08, the imaginary (Landau damping) term of *χ* for a kappa distribution is less than that of a Maxwellian, which the fitting algorithm compensates for by lowering the Maxwellian temperature. It is important to note that while the Maxwellian Landau damping term remains above the corresponding term for super-Gaussians at high *ξ* (corresponding to wavelengths far from the incident laser wavelength; see Fig. 1), the Maxwellian Landau damping term actually drops below the kappa Landau damping term at high *ξ*. This indicates that in a parameter regime where the VDF is a kappa distribution but the spectral peaks lie at high *ξ*, the error in the Maxwellian fitting could be reversed, which further illustrates the problems with applying the Maxwellian fitting algorithm to TS from non-Maxwellian plasmas.

Although the arbitrary fitting algorithm can perform better than the Maxwellian fitting algorithm in fitting TS spectra from non-Maxwellian VDFs, we have shown that the presence of non-Maxwellian VDF models can make the results of the fitting difficult to interpret without taking into account the details of how the different parameters of the non-Maxwellian VDF model correlate and affect the forward model. The use of the MCMC sampler on a synthetic distribution is one means by which we can study these correlations for given VDF models.

## VI. CONCLUSIONS

We have developed an arbitrary forward model that can compute integral-normalized TS spectra from arbitrary discretized electron and ion VDFs. The arbitrary forward model is benchmarked against a forward model that uses a tabulated plasma dispersion relation to compute the TS spectra for VDFs that are Maxwellian or linear combinations of Maxwellians. We discuss the numerical error associated with the arbitrary forward model and provide examples of how to minimize these effects.

The arbitrary forward model is used to implement an iterative fitting algorithm that accepts TS spectra as inputs and recovers plasma parameters defined by arbitrary user-defined VDF models. We show that the arbitrary fitting algorithm performs similarly to a Maxwellian TS fitting algorithm for Maxwellian VDFs but outperforms the Maxwellian fitting algorithm for non-Maxwellian VDFs. The arbitrary forward model can be used to fit the plasma parameters associated with any parameterizable non-Maxwellian VDF model that satisfies the basic TS assumptions of quasi-neutrality and an unmagnetized plasma (although this assumption could be relaxed with a suitable extension of the TS analytic model).

The fitting algorithm can be run with an MCMC sampler to estimate the PPD over the parameter space. This enables us to estimate the uncertainties of the fitted parameters and show that the numerical errors associated with the arbitrary forward model do not have a significant impact on the accuracy of the fitting. In the case of kappa VDFs, we also determine from the PPD that several parameters are linearly dependent (i.e., degenerate). The use of the arbitrary fitting algorithm to fit synthetic spectra from other non-Maxwellian VDF models could also reveal the correlated parameters in those models. Further work is needed to gain a better understanding of which VDF models contain these correlations and methods by which we can best mitigate them. Each VDF model could be examined separately using the tools we have developed in order to determine what parameter correlations are relevant in the PPDs. In addition, the computationally heavy calculation of *χ* and the resulting runtime increase could make the arbitrary forward model inconvenient to use in fitting experimental data. Further work should be performed in either optimizing the speed of the forward model or in designing other forward model schemes. For instance, forward models that are tailored for specific non-Maxwellian VDF models could be developed and benchmarked against our arbitrary forward model. Although these types of forward models would only be suited for those particular VDF models and thus inherently less generalized, it is possible that additional speed optimizations can be made in these cases.

## ACKNOWLEDGMENTS

We thank R. Follett for many valuable discussions related to this work. This work was supported by the U.S. Department of Energy (DOE) National Nuclear Security Administration (NNSA) under Award Nos. DE-NA0004033, DE-NA0003856, and DE-SC0020431, the University of Rochester, and the New York State Energy Research and Development Authority. This work was also supported by NASA under Grant No. 80NSSC19K0493.

This report was prepared as an account of work sponsored by an agency of the U.S. Government. Neither the U.S. Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof.

This research made use of PlasmaPy, version 2023.1.0, a community-developed open source Python package for plasma research and education (PlasmaPy Community *et al.*, 2023).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**B. C. Foo**: Formal analysis (lead); Investigation (lead); Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). **D. B. Schaeffer**: Conceptualization (equal); Funding acquisition (lead); Methodology (equal); Project administration (lead); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **P. V. Heuer**: Conceptualization (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*Plasma Scattering of Electromagnetic Radiation: Theory and Measurement Techniques*

*et al.*, (