A machine learning based Bayesian optimization solution to nonlinear responses in dusty plasmas

Machine learning (or deep learning) has recently become one of the hottest analysis techniques in the scientific world as application of this powerful numerical method has proven useful in solving problems across a wide range of fields. For example, convolutional neural networks (LeNet, AlexNet) Lecun et al. (1998); Krizhevsky et al. (2017) now fulfill object recognition tasks to a high degree of accuracy, recurrent neural networks (LSTM) Hochreiter and Schmidhuber (1997) are improving computational understanding of natural language, reinforcement learning agents are out-performing human experts in strategic decision making (AlphaGo) Silver et al. (2016, 2017), and Generative Adversarial Networks (GANs) Goodfellow et al. (2014) are showing the ability to create music, paintings and dialogue in a human manner. In addition to industrial applications, machine learning techniques are now also being applied to solve physics problems. Examples include the prediction of molecular atomization energies by employing regression models Rupp et al. (2012), the application of a neural network to solve quantum many-body problems Carleo and Troyer (2017), and crystallization recognition through the use of a shallow neural network Dietz et al. (2017).

In this paper, a machine learning-based method is applied to nonlinear response problems in dusty plasmas. Dusty plasmas Melzer et al. (1996); Fortov et al. (2005); Morfill and Ivlev (2009) are systems containing both weakly ionized gas and charged micron-sized dust particles. Due to the higher thermal velocities of electrons compared to ions, dust particles in a dusty plasma become negatively charged Goree (1994) in response to the frequent collisions between the plasma particles and the dust grain’s surface. Dust particle behavior in plasmas is determined by many factors, with the restoring confinement caused by the balance between the electrostatic force and gravity, the neutral gas drag, and particle-particle interactions between dust particles among the primary of these. Understanding the physics behind dust particle behavior (i.e., investigating these factors) is one of the most important tasks in dusty plasmas. One of the ways in which this can be fulfilled is by studying the response of the particles to external excitations Carstensen et al. (2011); Ding et al. (2019a). This is known as the nonlinear frequency response analysis, which has a wide application in mechanics, material science and nano science Peng et al. (2007); Cottone et al. (2009); Samanta et al. (2015).

Here, a Bayesian optimization framework Močkus (1975) is used to resolve a nonlinear response analysis Ivlev et al. (2000); Zafiu et al. (2001); Wang et al. (2002); Ding et al. (2019b) in a numerical manner in a dusty plasma. The undetermined coefficients in a dust grain’s equation of motion will be derived by optimizing the simulated motion to match that obtained from experimental results which will be compared to the analytic results from a multiple-scale perturbation method. Also, we will measure the nonlinearity to the 4th order in displacement, which helps correctly characterize the potential energy of particle in the plasma sheath but has not been investigated before .

It is important to note that this framework is not limited to nonlinear response analysis, but can also be applied to the more general case of physics problems where the experimental results can be reproduced by simulations. In these cases, undetermined physics quantities can be revealed efficiently (especially when the simulation is very computational expensive) by optimizing the simulations to experimental results in this Bayesian manner.

The experiment which will be discussed in this paper was conducted in a modified Gaseous Electronics Conference (GEC) RF reference cell (see Fig. 1) filled with argon gas. A single melamine formaldehyde (MF) particle having a diameter of 8.89 $\pm$ 0.09 $\mu m$ was inserted into a glass box (height: 20 mm, length: 18 mm, width: 18 mm) placed on the lower electrode which was powered at 13.56 MHz. The plasma power and pressure were fixed at 1.68 W and 40 mTorr, respectively. The MF particle was levitated in the plasma sheath region due to the balance between gravity and the electrostatic force produced by the negatively charged lower electrode. The dust particle was illuminated by a laser sheet (wavelength of 660 $nm$) with the resulting motion recorded by a high speed camera mounted at the side port of the cell at a rate of 500 fps. A primary amplitude-frequency response curve was measured by applying a sinusoidal excitation signal to the lower electrode with a fixed amplitude at various frequencies. Particle motion was recorded and then transformed into the frequency domain (FFT spectrum) using a Fourier transform at each value of the excitation frequency. The peak height of the FFT spectrum at the excitation frequency was measured, providing the primary response at this excitation. The secondary (super-harmonic) response to the excitation (a nonlinear response), can also be measured from the peak height of the FFT spectrum at twice the excitation frequency.

The motion of a single such particle levitating inside the plasma sheath under a vertical sinusoidal excitation can be modeled as a forced oscillator Ivlev et al. (2000),

$\displaystyle\ddot{x}+\mu\dot{x}+\omega^{2}x+\alpha x^{2}+\beta x^{3}=Fexp{(i\Omega t)}+c.c.,$

(1)

where $\mu$ is the neutral drag coefficient, $\omega$ is the restoring constant, $\alpha$ and $\beta$ are the second and third order derivatives of the restoring field, $\Omega$ is the frequency of the sinusoidal excitation, $F$ is the amplitude of the excitation (in units of acceleration) and c.c. stands for the complex conjugate. Usually, the effective restoring force experienced by the particle at equilibrium can be approximated as a linear function in displacement (i.e., $-\omega^{2}x$ where $\omega$ is considered the natural resonance frequency) under the assumption that the particle is levitating in a region where the sheath can be considered to exhibit a perfect parabolic sheath potential Tomme et al. (2000); Carstensen et al. (2012). Unfortunately, this approximation becomes invalid in most realistic situations, such as where charge fluctuations are considered Zafiu et al. (2001); Wang et al. (2002), or when the oscillation of the dust particle is so large that the sheath potential can no longer be considered parabolic. In this case, it is necessary to extend the restoring force to the nonlinear regime as $-\omega^{2}x-\alpha x^{2}-\beta x^{3}$ with terms higher than $O(x^{3})$ ignored for simplicity.

Based on Eq. 1, the particle motion $x(t)$ as a function of time under an excitation with frequency $\Omega$ can be simulated (or numerically solved) given a set of known parameters $\theta=\{\mu,\omega,\alpha,\beta,F\}$. In this case, the particle motion was simulated employing the velocity Verlet algorithm, which updates the position and velocity in each iteration as

$\displaystyle\begin{split}x(t+dt)=x(t)+v(t)dt+\frac{a(t)}{2}(dt)^{2},\\ v(t+dt)=v(t)+\frac{a(t+dt)+a(t)}{2}dt,\end{split}$

(2)

where $dt$ is the time step of the simulation, $v(t)$ is the velocity at time $t$ and $a(t)$ is the acceleration normalized by the particle mass at time $t$ as determined by Eq. 1

$\displaystyle a(t)=-\mu v(t)-\omega^{2}x(t)-\alpha x^{2}(t)-\beta x^{3}(t)+2Fcos(\Omega t).$

(3)

Following the same approach described above for measuring the amplitude-frequency response curves from experiment, the primary and secondary amplitude-frequency response curves can also be measured from the simulated particle motion $x(t)$ by varying the excitation frequency $\Omega$ over the range of excitain frequencies used in the experiment.

This allows a parameter set $\theta^{*}=\{\mu^{*},\omega^{*},a^{*},b^{*},F^{*}\}$ characterizing the properties of the dust motion (which depend on properties of the nearby plasma environment) to be determined by searching the parameter space for the optimal set of parameters that generates a simulated amplitude-frequency response curve resembling the experimentally measured amplitude-frequency response curve. This process can be quantified as

$\displaystyle\theta^{*}=\underset{\theta}{\operatorname{argmin}}(L(R_{e},R_{s}(\theta))),$

(4)

where $R_{e}$ represents the experimentally measured response curves and $R_{s}(\theta)$ represents the simulated response curves for a given set of parameters $\{\mu,\omega,\alpha,\beta,F\}$. $L(R_{e},R_{s})$ is a measure of the difference between the experimentally measured and simulated response curves. In order to quantify this difference, we define $L$ as a function $L:\theta=\{\mu,\omega,\alpha,\beta,F\}\mapsto\mathbb{R}$ that maps a set of parameters to a real value which measures the ‘distance’ between the experimentally measured and the simulated response curve as

$\displaystyle L(\theta)=\sum_{i=1}^{N}(\frac{r_{e}(\Omega_{i})-r_{s}(\Omega_{i},\theta)}{r_{e}(\Omega_{i})})^{2},$

(5)

where $r_{e}(\Omega_{i})$ and $r_{s}(\Omega_{i})$ are the experimentally measured and simulated response amplitudes at the excitation frequency $\Omega_{i}$, respectively. The summation is carried out over the excitation frequency span conducted in the experiment. The square operation is taken to ensure that $L(\theta)$ does not yield negative values, which guarantees the existence of minimal points. It is important to mention that this type of loss function $L(\theta)$ may not be well-defined everywhere. For unrealistic parameters sets, i.e., sets that either have no physical meaning or are not suitable for describing the condition of the plasma sheath, this simulation of nonlinear response curves diverges, resulting in an undefined distance function. In these cases, a large value is assigned to the distance function (e.g., $L=10^{5}$) in order to ensure optimization success.

As can be seen from Eq. 5, calculation of the loss function $L(\theta)$ for even one set of parameters requires multiple simulations of the particle’s motion, i.e., one for each excitation frequency used for measuring the response curve in the experiment. For example, the response curve shown in Fig. 2 requires 71 independent simulations to calculate the distance function for just one set of parameters. This is computationally expensive and, as such, a minimization of the distance function $L(\theta)$ based on a random search of the parameter space $\theta$ is infeasible.

Therefore, this loss function is minimized employing a Bayesian optimization. This technique has shown great promise in machine learning, especially for the fine tuning of neural networks for model selection. Now, instead of randomly searching the parameter space and then conducting simulations for each set, only those parameter sets selected in a Bayesian manner are simulated. A surrogate function $f$ is introduced to model the distribution of the value of the loss function $L(\theta)$. The posterior distribution of this surrogate function at the parameter $\theta$ given the data observed $\mathcal{D}_{1:t}=\{\theta_{1:t},L(\theta_{1:t})\}$ can now be derived using the Bayes law

$\displaystyle p(f|\theta;\mathcal{D}_{1:t})=\frac{p(\theta|f;\mathcal{D}_{1:t})p(f;\mathcal{D}_{1:t})}{p(\theta;\mathcal{D}_{1:t})},$

(6)

where Tree-structured Parzen density estimators Bergstra et al. (2011) (a generative model) are used to model the likelihood function $p(\theta|f;\mathcal{D}_{1:t})$ defined as

$$p(\theta|f;\mathcal{D}_{1:t})=\begin{cases}l(\theta),\text{if }f<f^{*}\\ g(\theta),\text{if }f\geq f^{*}.\end{cases}$$

(7)

In this likelihood function, $l(\theta)$ and $g(\theta)$ are non-parametric Parzen density estimators (i.e., Gaussian mixtures) of the likelihood for the data to have a value of the loss function $L(\theta)$ smaller and greater than a threshold $f^{*}$, respectively. As such, the marginal distribution of the parameter set given the observed data set $\mathcal{D}_{1:t}$ (the denominator of Eq. 6) can in turn be calculated as

$\displaystyle\begin{split}p(\theta;\mathcal{D}_{1:t})=&\int_{-\infty}^{\infty}p(\theta|f;\mathcal{D}_{1:t})p(f;\mathcal{D}_{1:t})df\\ =&(l(\theta)-g(\theta))\int_{-\infty}^{f^{*}}p(f;\mathcal{D}_{1:t})df+g(\theta).\end{split}$

(8)

The criteria for exploring the overall parameter space is to choose the next set of simulation parameters that maximizes the expected improvement $\mathbb{E}[max(f^{*}-f,0)]$ Jones et al. (1998) as

$\displaystyle\begin{split}\theta_{t+1}=&\underset{\theta}{\operatorname{argmax}}\int_{-\infty}^{\infty}\operatorname{max}(f^{*}-f,0)p(f|\theta;\mathcal{D}_{1:t})df\\ =&\underset{\theta}{\operatorname{argmax}}\frac{\int_{-\infty}^{f^{*}}(f^{*}-f)p(f;\mathcal{D}_{1:t})df}{\frac{g(\theta)}{l(\theta)}(1-\int_{-\infty}^{f^{*}}p(f;\mathcal{D}_{1:t})df)+\int_{-\infty}^{f^{*}}p(f;\mathcal{D}_{1:t})df}\\ =&\underset{\theta}{\operatorname{argmax}}\frac{l(\theta)}{g(\theta)},\end{split}$

(9)

where the last equation holds since the cumulative distribution $\int_{-\infty}^{f^{*}}p(f;\mathcal{D}_{1:t})df$ is strictly less than 1. As shown, this result is not affected by the exact form of the prior $p(f;\mathcal{D}_{1:t})$. As such, the next parameter set whose distance function will be simulated is chosen to maximize the quotient of the Parzen density estimators $l(\theta)/g(\theta)$. As each new simulation is conducted, the data set $\mathcal{D}$ will be updated with the new simulated data points (e.g., the data set $\mathcal{D}_{1:t}$ is updated to $\mathcal{D}_{1:t+1}$ by adding a new simulated point $\{\theta_{t+1},L(\theta_{t+1})\}$). The posterior distribution of the surrogate function Eq. 6 will also be updated accordingly, eventually resembling the behavior of the real distance function Eq. 5. Notice that the threshold $f^{*}$ in Eqs. 7 and 9 are the up-to-date optima (i.e., the minima of the loss function updated to the data $\mathcal{D}$).

Since secondary responses (as nonlinear responses) are very sensitive to nonlinear terms, i.e., $\alpha x^{2}$ and $\beta x^{3}$, while the primary responses are more sensitive to the linear terms, it is necessary to minimize the loss functions for both primary and secondary responses simultaneously. One simple way of achieving this is to minimize a weighted sum of these two loss functions rather than minimizing them individually (e.g., $L=L_{p}+0.05L_{s}$).

Figure. 2 shows the Bayesian-optimized simulated primary response curves (dashed red curves) of a single dust particle levitated in the plasma sheath in the GEC RF reference cell at a plasma power of 1.68 Watts and pressure of 40 mTorr. The corresponding secondary response curves are shown in the subplots. Particles excited under excitation amplitudes of 1.0 V and 1.5 V are plotted in Fig. 2a and Fig. 2b, respectively. As shown, the optimized response curves (dashed red curves simulated according to Eq. 1) resemble the experimentally measured responses curves (solid black curves) in both the primary and secondary regions. Also, note that the spring softening phenomenon (i.e., the nonlinear phenomenon that results in the primary resonance peak being ‘bent’ in the low frequency direction) becomes more obvious as the excitation amplitude increases (Fig. 2a).

The corresponding optimized parameters obtained from model 1 (Eq. 1) are calculated as the average over five independent trials of the optimizing experiment and their values are listed in Table 1 (method ‘Model 1’), with the corresponding Coefficient of Variations (CV) shown in parentheses. Notice that the sign of the coefficient of the quadratic nonlinearity $\alpha$ is irrelvant in this response analysis since it only changes the direction of the asymmetric motion of the dust particle. Even though there is a large variation in randomness in this Bayesian search of the parameter space, the optimizing experiment converges to yield consistent results as evidenced by the low CV. The relatively high CV for the parameter $\beta$ (the coefficient of the cubic nonlinearity) is due to the fact that the response curves are robust in responding to variation of nonlinearities of higher order. As such, a small variation in $\beta$ will not significantly perturb the entire response curve.

Fig. 3 shows the loss (Eq. 5) as a function of the number of iterations for the dust particle excited at both 1.0 V (a) and 1.5 V (b), each of which has five independent experimental trials. As shown, the loss values decrease rapidly after the first several iterations, reaching convergence after a few hundred iterations. (This again indicates the efficiency of the presented Bayesian optimization method in exploitation of the parameter space.) However, in order to boost overall accuracy and ensure wider exploration of the parameter space, we conducted a large number of iterations. The observed higher convergency loss value for the dust particle under 1.5 V excitation ($\approx 4.3\times 10^{-3}$) as compared to that under 1.0 V excitation ($\approx 2.2\times 10^{-3}$) can be attributed to the increased difficulty of capturing the spring softening phenomenon as the excitation amplitude becomes larger. (Compare Fig. 2a and Fig. 2b.)

The parameters can also be derived analytically by solving the equation of motion (Eq. 1) employing the multiple-scale perturbation method. The details of this method are given in refs Nayfeh and Mook (1979); Ding et al. (2019b), with the main results needed for the analysis given below.

Assuming an external excitation at a frequency of approximately half that of the oscillator resonance frequency $\omega$, i.e., $\Omega\approx\frac{1}{2}\omega$, the solution to Eq. 1, to first order, yields

$\displaystyle\begin{split}x(t)=&\frac{F}{\omega^{2}-\Omega^{2}}cos(\Omega t)\\ &-\frac{\alpha F^{2}}{4\omega(\omega^{2}-\Omega^{2})^{2}\sqrt{\frac{\mu^{2}}{4}+(2\Omega-\omega)^{2}}}sin(2\Omega t-\phi),\end{split}$

(10)

where $\phi$ is the shifted phase which is dependent on the excitation frequency as $\phi=arctan(\frac{4\Omega-2\omega}{\mu})$. The parameters $\omega$, $\mu$ and $\alpha$ are determined using the experimentally measured secondary response curve fitted to the steady state theoretical secondary response $\alpha F^{2}/4\omega(\omega^{2}-\Omega^{2})^{2}\sqrt{\frac{\mu^{2}}{4}+(2\Omega-\omega)^{2}}$. Considering an external excitation having a frequency approximately equal to the oscillator resonant frequency, i.e., $\Omega\approx\omega$, the solution to Eq. 1, to first order of approximation yields

$\displaystyle\begin{split}x(t)=A(\Omega)cos(\Omega t-\phi^{\prime}),\end{split}$

(11)

where the shifted phase is $\phi^{\prime}=\Omega-\omega-\beta$. By eliminating the secular term appearing in the equation of motion to second order of approximation, the steady state theoretical primary response $A(\Omega)$ can now be derived as

$\displaystyle\begin{split}\frac{F^{2}}{4\omega^{2}}=(\frac{A\mu}{2})^{2}+[(\frac{9\beta\omega^{2}-10\alpha^{2}}{24\omega^{3}})A^{3}-(\Omega-\omega)A]^{2}.\end{split}$

(12)

By fitting the corresponding experimentally measured primary response curve to Eq. 12, the parameter space for the cubic nonlinearities $\beta$ and $F$ can now be investigated. The parameters obtained in this way are shown in Table 1 (method ‘Multiple-scale’). As shown, the parameters measured from the Bayesian optimization are consistent with those measured from the multiple-scale perturbation, except for the value of $\beta$ under 1.5 V excitation (with approximately 57.6$\%$ percent difference). Notice that the measurements of the parameters from the multiple-scale perturbation serve as a benchmark and should not be considered as true values since they are derived and are precise only to the first order of approximation.

Due to limitations of the multiple-scale perturbation method, any extension of the model (Eq. 1) (e.g., including higher order nonlinearities of either displacement $x$ or velocity $\dot{x}$) and derivation of the corresponding approximate solutions would be tediously complicated. However, the Bayesian optimization scheme described here allows this process to be simplified greatly. As an example, the model is extended to include an additional nonlinearity of higher order in displacement $x$,

$\displaystyle\ddot{x}+\mu\dot{x}+\omega^{2}x+\alpha x^{2}+\beta x^{3}+\gamma x^{4}=Fexp{(i\Omega t)}+c.c..$

(13)

Applying the Bayesian method, the optimized parameters are measured, with the results shown in Table 1 (with method ‘Model 2’) and the corresponding simulated response curves shown in Fig. 2 (dashed blue curves). Again, the corresponding CV for the five independent trials are shown in parentheses.

By considering nonlinearities to fourth order, the primary response curves (dashed blue line) in Fig. 2a more closely resemble the spring softening behavior than do those considering nonlinearities to third order (dashed red line). Also, the loss is further reduced, reaching $1.6\times 10^{-3}$ for a 1.0-V excitation and $2.3\times 10^{-3}$ for a 1.5-V excitation, as compared to values based on the model with third order nonlinearities (Eq. 1) which results in loss values of $2.2\times 10^{-3}$ and $4.3\times 10^{-3}$, respectively. This indicates a closer match of the simulated response curves to the experimentally measured ones. After introducing nonlinearities to the fourth order, the measured drag coefficient $\mu$, excitation amplitude $F$, and the coefficient of the quadratic nonlinearity $\alpha$ more closely approach the values measured from the multiple-scale perturbation. However, it is noted that the coefficient for the cubic nonlinearities, $\beta$ exhibits a large deviation. Considering the conditions for the existence of the spring softening effect (seen from Eq. 12)

$\displaystyle\begin{split}9\beta\omega^{2}-10\alpha^{2}<0,\end{split}$

(14)

the critical value of $\beta$ for the existence of the spring softening phenomenon can be derived as $\beta<\beta_{c}\approx 8.1\times 10^{-4}$ (by taking into consideration the fact that the measured values for the coefficient of the quadratic nonlinearity $\alpha$ are consistent in both models given by Eq. 1 and Eq. 13). In this case, a large value of $\beta$ as measured for both a 1.0-V excitation and a 1.5-V excitation, based on the model represented by Eq. 13, seems to violate the condition of the existence of the spring softening phenomenon (Eq. 14). However, Eq. 14 is derived from only the first order of approximation in the multiple-scale perturbation. Thus, although a response curve simulated with a parameter $\beta$ violating the condition given by Eq. 14 still reveals the spring softening phenomenon, this indicates a limited ability and accuracy of the multiple-scale perturbation method to explain nonlinear responses since it ignores higher order nonlinearities. In order to accurately determine the coefficient of the cubic nonlinearities $\beta$ (an important factor characterizing the nonlinearity of the plasma sheath Ivlev et al. (2000); Zafiu et al. (2001)), the effects from higher order nonlinearities are important and should not be ignored. Fig. 4a shows the effective restoring potential energy $\Phi$ of the particle (divided by the particle mass) in the vicinity of its equilibrium position for both Model 1 (red) and Model 2 (blue). The difference in restoring potential energy is observable for these two models.

With the coefficient for the fourth order nonlinearity available, the change in the electric field and the grain charge at varying levitation positions can be further investigated. By considering an expansion in the electric field $E$ and grain charge $Q$

$\displaystyle\begin{split}&E=E_{0}+E_{1}x+E_{2}x^{2}+E_{3}x^{3},\\ &Q=Q_{0}+Q_{1}x+Q_{2}x^{2}+Q_{3}x^{3},\end{split}$

(15)

the electrostatic force can be written as

$\displaystyle\begin{split}F_{stat}=&(E_{0}Q_{0})+(E_{0}Q_{1}+E_{1}Q_{0})x\\ &+(E_{0}Q_{2}+E_{1}Q_{1}+E_{2}Q_{0})x^{2}\\ &+(E_{0}Q_{3}+E_{1}Q2+E_{2}Q_{1}+E_{3}Q_{0})x^{3}\\ &+(E_{1}Q_{3}+E_{2}Q_{2}+E_{3}Q_{1})x^{4}.\end{split}$

(16)

The coefficents of the polynomial in this expansion are related to the corresponding coefficents in Eq. 13. By assuming a linear electric field (i.e., $E_{2}=E_{3}=0$), the nonlinearities in charge $Q$ can be explored to the third order ($Q_{3}$) with $\gamma$ provided by the Bayesian optimzation approach (Model 2), while without $\gamma$, the nonlinearities in charge can only be explored up to the second order ($Q_{2}$). Fig. 4b shows the grain charge in the vicinity of the equilibrium position for the third order polymonial ($Q=Q_{0}+Q_{1}x+Q_{2}x^{2}+Q_{3}x^{3}$) and second order polymonial ($Q=Q_{0}+Q_{1}x+Q_{2}x^{2}$) expansion in red and blue, respectively. As a reference, a linear charge model ($Q=Q_{0}+Q_{1}x$) is also shown in black. The corresponding linear electric field (divided by the particle mass) for each charge expansion are shown in Fig. 4c. The equilibrium charge $Q_{0}\approx 1.4\times 10^{4}e$ was estimated by using the levitation position comparison method Ding et al. (2016). In this method, a vertically aligned two particle pair is formed, and the difference in the levitation position for the upstream particle with and without the presence of the downstream particle (where the downstream particle is knocked out of the system using a laser pulse), is measured. As shown, the third order polymonial charge model predicts a weaker charge reduction in the downstream region, but a stronger charge accumulation in the region above the equilibrium position.

Beyond assuming a linear electric field, we can also investigate the nonlinear expansions for the E-field and the grain charge simultaneously (i.e., $E=E_{0}+E_{1}x+E_{2}x^{2}$ and $Q=Q_{0}+Q_{1}x+Q_{2}x^{2}$). However, due to a lack of constraints, this investigation can only be explored to the second order for both electric field and grain charge even with $\gamma$ known. Fig. 4b and Fig. 4c show the result of this charge model and the corresponding nonlinear electric field as dashed lines. As shown, both the grain charge and the electric field can be very different from the cases where the electric field is assumed to be linear in the downstream region. This indicates a resonable assumption of a linear electric field in a close vicinity of the equilibrium position. It is instructive to compare these results against the usual linear models for both the particle charge and electric field. With the assumption of a linear electric field, the charge varies considerably from the linear charge model in the upstream direction. Conversely, with the assumption of a quadratic electric field, it is seen that the charge varies significantly from the linear charge model in the downstream direction. Future experiments may be designed to determine which of these models is correct.

In conclusion, a nonlinear response analysis for dust particles in plasma was provided employing a machine learning based method. An efficient technique for optimizing the comparison between numerically simulated and experimentally measured response curves by searching the parameter space in a Bayesian manner was described. Using this approach, the physical parameters characterizing the plasma conditions can be derived. The nonlinearity of the response was determined to the fourth order, which is necessary in order to accurately determine the the coefficients for lower-order nonlinearities, as well as to correctly characterize the potential energy of the particle in the sheath. Beyond the field of dusty plasmas, the proposed framework provides a general method for measuring physical quantities by optimizing simulation parameters to match experimental observations in an efficient manner, especially when the simulation is computationally expensive.

Support from NASA grant number 1571701, NSF grant number 1740203 and NSF grant number 1707215 is gratefully acknowledged.