Deriving Heating Model from E Field Fluctuations

Previously I have tried to consider the Hockney calculations regarding the E-Field Fluctuations, but I did not fully examine how this would turn into a heating model. Here, I will finish the calculation. In the Hockney 1971 paper, they state a formula for the time-averaged squared electric field fluctuations (which comes from Plasma Kinetic Theory by Montgomery and Tidman Ch. 8.2)

\begin{equation} \frac{\langle E^2 \rangle}{8 \pi} = \frac{k_B T}{2} \int \int_{-\infty}^{+\infty} \frac{d k_x d k_y}{(2 \pi)^2} \frac{1}{(1 + (k_x^2 + k_y^2) \lambda_D^2)} \end{equation}

which is in units of energy per unit volume. We can convert this integral to polar coordinates by letting $k = \sqrt{k_x^2 + k_y^2}$ and integrating with respect to $k \cdot dk \cdot d\phi$. Also, we can convert the expression from CGS to SI Units

\begin{equation} \frac{\epsilon_0 \langle E^2 \rangle}{2} = \frac{k_B T}{4 \pi} \int_0^\infty \frac{k dk}{1 + k^2 \lambda_D^2} \end{equation}

To solve this integral, we let $u \equiv k \lambda_D$ to get

\begin{equation} \langle E^2 \rangle = \frac{k_B T}{2 \pi \epsilon_0 \lambda_D^2} \int_0^\infty \frac{u du}{1 + u^2} \end{equation}

which diverges at infinity, so we set limits to the integral 0 and $u_{max}$ which correspond to the maximum relevant wavenumber in our simulation. This can be solved analytically as

\begin{equation} \langle E^2 \rangle = \frac{k_B T}{4 \pi \epsilon_0 \lambda_D^2} \text{Log}(1 + u_{max^2}) \label{eq:e2} \end{equation}

Here, $u_{max} = 2 \pi \frac{\lambda_D}{\Delta x} R$ where $R$ is some multiplicative factor that tells us where to implement our cutoff. $R = 1$ corresponds to saying that the maximum wavenumber considered is limited the the cell size: $k_{max} = \frac{2 \pi}{\Delta x}$.

Also, from the Hockney 1971 Paper, they give an estimate of how the squared electric field fluctuations contribute to the kinetic energy increase of one macro particle

\begin{equation} \Delta KE = \frac{1}{2} m \langle \Delta v^2 \rangle = \frac{q^2}{m} \langle E^2 \rangle \tau_{corr} t \end{equation}

where we can identify $\tau_{corr} = \frac{2 \pi}{\omega_{pe}}$ with the plasma period and we already have $\langle E^2 \rangle$. Additionally, we can express the charge of one macroparticle as $q = e \frac{N}{N_{mac}} = \frac{e n}{n_{mac}} = \frac{e n \Delta x^2}{n_{ppc}}$ Combining this all together, we will get

\begin{equation} \Delta KE = (\frac{e}{m_e}) \frac{e n \Delta x^2}{n_{ppc}} \frac{k_B T}{4 \pi \epsilon_0 \lambda_D^2} \text{Log}(1 + u_{max}^2) \frac{2 \pi}{\omega_{pe}} t \end{equation}

now, we can substitute the expressions for the plasma oscillation frequency and debye length to arrive at

\begin{equation} \Delta KE = \frac{n_e^{3/2} \Delta x^2}{n_{ppc}} \text{Log}(1 + u_{max}^2) \sqrt{\frac{e^6}{4 m_e \epsilon_0^3}} \end{equation}

From here, I can do the same procedure I did with the arber formula and let $n \rightarrow 10^{29} n_{23}$, $\Delta x \rightarrow 10^{-9} \Delta x_{nm}$, and $t \rightarrow 10^{-12} t_{ps}$ to get

\begin{equation} \frac{\Delta T_{eV}}{\Delta t_{ps}} = (1.07 \times 10^7 eV) \frac{n_{23}^{3/2} \Delta x_{nm}^2}{n_{ppc}} \text{Log}(1 + u_{max}^2) \end{equation}

which, to me, does not line up with the heating constant we see in the Arber formula, so I will just call the constant out front $C_H$. Additionally, I will add another hyper-parameter that multiplies the debye ratio called $R$ and when $R=1$, we would arrive at the final formula:

\begin{equation} \frac{\Delta T_{eV}}{\Delta t_{ps}} = C_H \frac{n_{23}^{3/2} \Delta x_{nm}^2}{n_{ppc}} \text{Log}(1 + 4 \pi^2 (R \frac{\lambda_D}{\Delta x})^2) \end{equation}

which is identical to the Arber-Model, with the exception of the Log Term. When the cell size resolves the debye length ($\lambda_D > \Delta x$), the log term is very slowly growing as $\Delta x$ gets smaller, but this is overshadowed by the $\Delta x^2$ term outside of the Log. As a result, this would not behave much differently than the Arber formula when the cell size is small, which is in the regime where there is not much numerical heating anyway.

When the Debye Length is very poorly under-resolved ($\lambda_D \ll \Delta x$), we can simplify the expression as

\begin{equation} \frac{\Delta T_{eV}}{\Delta t_{ps}} = C_H \frac{n_{23}^{3/2} \Delta x_{nm}^2}{n_{ppc}} 4 \pi^2 (R \frac{\lambda_D}{\Delta x})^2 = 4 \pi^2 R^2 C_H \frac{n_{23}^{3/2}}{n_{ppc}} \lambda_{D, nm}^2 \end{equation}

which is now no longer dependent on the cell size $\Delta x$. This is consistent with our observation that after a certain cell size is reached, the amount of heating becomes relatively constant.

Comparison with My Simulations

I restested this model (as well as the alpha and ratio models) with my data that I generated described here. Just a slight notational difference: now I call the heating constant out front $C_H$ which does not include the factor of $n_{ppc} = 25$ anymore and I call the multiplier in both the ratio and hockney models $R$.

First, I looked at the Arber Model which had a MSE of 20.2. For a more fair comparison, If I allow a hyper-parameter of a heating constant $C_H$ (which is 2 in Arber’s Paper for the current settings) I get a slightly better fit with $C_H = 1.8$ which reduced the MSE to 19.6 (improvement is minimal, which shows Arber’s Heating Constant to be generally consistent with my data).

Next, I repeated the Alpha and Ratio Model (here I gave similar findings) which had MSE of $6.43$ and $3.69$ respectively. (Before I restricted the Alpha Model to having $\alpha$ being a half-integer, so the estimate is slightly better now that we have $\alpha = 0.66$).

Finally, I fitted my new model based on E-Field Fluctuations in the Hockney Paper, which I will call the Hockney Model, that had a $MSE = 3.25$ which is similar to the ratio model as expected. To enable a more direct comparison with the Arber Model, I set $R = 1$ in the Hockney Model so that there is only one hyper-parameter of the heating constant and got a MSE = 12.64 which is definitely an improvement over the Arber-Model’s 19.6.

Below, I summarize these results a table. I have two sets of model fits: the first corresponds to minimizing the MSE with absolute errors and the second considers MSE with percentage error.

Next, I made some plots comparing the amount of heating from each model vs the amount of heating in the simulations I ran. I made plots for each sim at a constant debye ratio (13 sims per debye ratio, 7 different debye ratios). Here, I am just considering the parameters for the models using the normal MSE with absolute error.

Just to get a feel for how much heating is expected, I plotted the prediction for the models compared to the actual amount of heating in the models.

We notice that the amount of heating definitely varies on the simulation within a constant debye ratio. This shows that the debye ratio is not enough to predict the amount of numerical heating. Additionally, simulation index 8 is the one that has the highest temperature and temperature ($900 eV$, $n_{23} = 9$) which makes sense why it would have a higher temperature gain than the others. Additionally, we definitely see that the Arber Model is inadequate.

To get a better look at the model comparisons, I also plotted these deviations as a percentage difference of the simulation result.

Here, it is more clear that the ratio and hockney model perform similarly and that they do a decent job of predicting the amount of numerical heating in a simulation. However, we still see a trend in that the model under-predicts for larger debye ratio and slightly over-predicts for low debye ratio.

To see how changing the ratio multiplier $R$ affects the model, I calculated the MSE and $C_H$ for a given $R$ across a range of $R$ from $0 \rightarrow 20$. Below I have the MSE vs R

And here, I do the same graph, but with the heating constant $C_H$ on the x axis

From this parameter scan, I get that the optimal parameters are $C_H = 7.4$ and $R = 3.7$ which has a $MSE = 3.18$ (this is slightly different than what I obtained before just because my parameter scan this time is much tighter). Since the B-Spline 3 Shape Function has an effective width of $4 \Delta x$, I wonder if this could be related.

Comparison with Ricky’s Simulations

Below, I do the same comparisons with Ricky’s No Laser Simulation.

Ricky’s simulations are much harder to parse. There are fewer sims overall, and the sims that have a larger number of particles per cell cannot be used because they don’t heat up over the course of the simulation. However, it is clear that the Arber Formula does not work well at all for this set. And, at the very least, the Hockney 1 Parameter Model performs better than the Arber 1 Parameter Model (just heating constants out front) which indicates that the Hockney Model is still an improvement over the Arber Model.

I can do a paramter scan of the multiplier in the Hockney Model to determine an optimal multiplier

Here, the mean RMS error fluctuates more which would mean that the model is very sensitive to both parameters we pick. As a result, the model doesn’t appear to be approximating the data well, but instead just overfitting two hyper-parameters.

Things to Do

• Maybe Try redoing sims for triangle and tophat shapefunctions