Electric Field Fluctuations

From the 1971 Hockney Paper, there was also a section that described the time average of the squared E-field fluctuations with respect to the particle kinetic energy. This is important because we know that the squared field fluctuations are related to numerical heating.

Hockney’s calculation starts from the work by Montgomery and Tidman Ch. 8. I borrowed this book from the library and tried to understand where he was coming from, but it didn’t really make any sense to me. It involved auto-correlation functions and moving to the frequency domain among other things.

Hockney arrives at

\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}

and in 2D, $k_B T = mv_{th}^2$ and if we let $u = k \lambda_D$ (where $k$ is the wavenumber and not the boltzmann constant) we can convert this to an integral in polar coordinates and arrive at

\begin{equation} \frac{\langle E^2 \rangle / 8 \pi}{n m v_{th}^2} = \frac{1}{16 \pi} \text{Log} (\frac{1 + u_{max}^2}{1+u_{min}^2}) \frac{1}{n \lambda_D^2} \end{equation}

In our PIC sims, $n$ is the macro-particle density per unit area and $u_{max} = \pi (\lambda_D / \Delta x)$ and $u_{min} = (2 \pi / N_x)(\lambda_D / \Delta x)$ would be the limits of the integral based on the relevant length scales in the simulation grid. If we plug in these limits, we would get

\begin{equation} \frac{\langle E^2 \rangle / 8 \pi}{n m v_{th}^2} = \frac{1}{16 \pi} \text{Log} (\frac{1 + (\pi r)^2}{1+(2 \pi r / N_x)^2}) \frac{1}{n \lambda_D^2} \end{equation}

where $r = \lambda_D / \Delta x$ is the debye ratio.

However, I have an issue with this equation. Let’s say we fix the Debye Length and number density to some value. Then, we can plot these E-field fluctuations (ignoring the constants out front) as a function of only the debye ratio.

Above, I have plotted this function where the variable $x$ is the debye ratio. We can see as the debye ratio increases, the amount of heating also increases monotonically. This is not what we observe however in our sims, so I don’t think that this would be a good model. However, we said $n$ is the density per unit area, but if we use that as the density per area of one grid cell, then we would get a plot like

Which has similar behavior to what we expect. When the debye ratio is close to zero, the heating plateaus to some amount instead of increasing to infinity. And the heating goes as 1 over the square of the debye ratio when the ratio is close to 1. Also, instead of the Log formula, we can use $(x+1)^{-2}$ where x is proportional to the debye ratio. Below, I have plotted these two as a rough comparison:

The reason I compare thse two formulas is that they have the same functional form for $x \gg 1$ and $x \ll 1$. So, if we were to write out the full formula for the heating for both of these sims, we would get

\begin{equation} \frac{dT}{dt} \propto \text{Log} (\frac{1 + (\pi r)^2}{1+(2 \pi r / N_x)^2}) \frac{1}{n_{ppc} r^2} \end{equation}

and

\begin{equation} \frac{dT}{dt} \propto (\frac{\Delta x}{\Delta x + \lambda_D} )^2 \frac{1}{n_{ppc}}= (\frac{1}{1 + r} )^2 \frac{1}{n_{ppc}} \end{equation}

where we would just need to fit some specific constant of proportionality.

Model Comparisons

After implementing these models, I realized that the models as is would definitely not work. That is because both of these new models are only dependent on the debye ratio and nothing else. I would group my simulations according to constant debye ratio, so these models would just naively predict a constant value for the heating for some debye ratio. However, we know this type of model is only adequate when the debye ratio is really small. When the debye ratio is closer to 1, we know this constant heating is not what is going on.

For example, I have the absolute amount of heating in (eV) displayed below for the four models (when $\Delta x = 3$ nm)

This was the set of simulations that gave the best results for the $\alpha$-model and we can clearly see that the log or ratio model is definitely not adequate.

Modified Arber Formula

Previously I showed where the Arber Formula came from where the heating time was given by the Hockney 1971 paper as

\begin{equation} \tau_H \approx 200 \pi \frac{n_{ppc}}{\omega_{pe}} ( \frac{\lambda_D}{\Delta x} )^2 \end{equation}

I did figure out how to correctly convert the units to the form Arber has in eq. 30 in his paper. First, solve for the rate of change of temperature and substitute the expressions for the debye length and plasma frequency

\begin{align} \frac{dT}{dt} &= \frac{T_0}{\tau_H} \\ &= T_0 \sqrt{\frac{n e^2}{m \epsilon_0^3}} \frac{1}{n_{ppc}} \frac{\Delta x^2}{200 \pi} \frac{n e^2}{\epsilon_0 k_B T_0} \\ &= \sqrt{\frac{e^6}{m \epsilon_0^3}} \frac{1}{200 \pi k_B} \frac{n^{3/2} \Delta x^2}{n_{ppc}} \\ &= 1.879 \times 10^{-5} \frac{n^{3/2} \Delta x^2}{n_{ppc}} \end{align}

Then we make the following substitutions:

• $\Delta x \rightarrow \Delta x_{nm} \times 10^{-9}$
• $n \rightarrow n_{23} \times 10^{29}$ (power of 29 is b/c $n_{23}$ is in reference to $cm^{-3}$ and not $m^{-3}$)

And this becomes

\begin{align} \frac{dT}{dt} &= 1.879 \times 10^{-5} 10^{29 \times 1.5} 10^{-9 \times 2} \frac{n_{23}^{3/2} \Delta x_{nm}^2}{n_{ppc}} \\ &= 5.942 \times 10^{20} \frac{n_{23}^{3/2} \Delta x_{nm}^2}{n_{ppc}} \end{align}

which has units of kelvin per second. If we convert this to electron volts per picosecond, we would get

\begin{align} \frac{dT_{eV}}{dt_{ps}} &= 5.942 \times 10^{20} \frac{n_{23}^{3/2} \Delta x_{nm}^2}{n_{ppc}} \frac{K}{s} \times \frac{10^{-12} s}{1 ps} \times \frac{1 eV}{11605 K} \\ &= 51201 \frac{n_{23}^{3/2} \Delta x_{nm}^2}{n_{ppc}} \end{align}

which lines up with his listed heating coefficient of $\alpha_H \sim 5 \cdot 10^4$. If instead I used the following for the heating time

\begin{equation} \tau_H \approx 200 \pi \frac{n_{ppc}}{\omega_{pe}} ( \frac{\lambda_D + \Delta x}{\Delta x} )^2 \end{equation}

I can show that this reduces to an Arber-Like formula of

\begin{equation} \frac{dT_{eV}}{dt_{ps}} \approx 30 \frac{T_0 \sqrt{n_{23}}}{n_{ppc}} \frac{1}{(\lambda_D / \Delta x + 1)^2} \end{equation}

where now, the temperature dependence $T_0$ does not cancel out. Also, the debye length is a function of density $n$ and temperature $T_0$, but I chose to leave it in this form so that the equation would not look too complicated. If we want to make a direct comparison to the Arber formula, we can use

\begin{equation} \frac{dT_{eV}}{dt_{ps}} = \frac{\alpha_H}{8} \frac{T_0 \sqrt{n_{23}}}{n_{ppc}} \frac{1}{(\lambda_D / \Delta x + 1)^2} \end{equation}

which has the same approximate numerical prefactor as the Arber Formula if the +1 in the denomenator were to go away (actually not sure about this one).

Characteristic Quantities of Simulation

Using the Heating Time Formula Above, we can compute an estimate of the scale of the heating time $\tau_H$.

\begin{equation} \tau_H = 200 \pi \frac{n_{ppc}}{\omega_{pe}} \frac{5 \times 10^4}{\alpha_H} \end{equation}

which would be 0.9 ps for the Hockney Paper and 22 ns for EPOCH with spline and current smoothing. This tells us that the amount of heating when the debye length is on the order of the grid spacing is neglible for the EPOCH sims. This also tells us that the ratio model would routinely under-predict the amount of heating (compared to the arber formula) unless we give it a heating constant that would be greater than the arber formula.

Additionally, if we are considering a simulation with density $n_{23} = 1$ and $\Delta x_{nm} = 3$, we would find the plasma frequency as

\begin{equation} \omega_{pe} = 1.78 \times 10^{16} \text{1/s} \end{equation}

which corresponds to a value of $\omega_{pe} \Delta t \sim 0.18$. We know there is an instability when $\omega_{pe} \Delta t > 2$, so this may only be an issue for some of the larger $\Delta x$ and $n_{23}$ simulations like when $\Delta x_{nm} > \sim 12$ nm.

Things to Do

• Run Ratio Model for Different Grid Sizes
• See what Modifications to Ratio Model would make a better model
• Think of theoretical basis for heating cutoff at certain debye ratio.
• Ask Ricky for his sims without the target to see if I can extract any info about that.