Electrostatic linear plasma dispersion relation
$$
1 + \Sigma_s \frac{\Pi_s^2}{\omega}\int_{-\infty}^{\infty} \frac{\partial f_s(v_\parallel)}{\partial v_\parallel}\frac{v_\parallel}{\omega-k_\parallel v_\parallel} dv_\parallel=0
$$
where, in SI units:
-
$\Pi_s=\sqrt{\frac{q^2 n}{m\epsilon_0}}$ is the plasma frequency of species $s$,
-
$q$ is the charge of the physical particle (not computational macro-particle) of species $s$,
-
$m$ is the mass the physical particle
-
$n$ is the number density of the species,
-
$\epsilon_0$ is the electric permittivity,
-
$\omega$ is the complex frequency
-
$v_\parallel$ is the one-dimensional velocity coordinate
-
$f_s(v_\parallel)$ is the distribution function of species $s$ and is normalised such that $\int_{-\infty}^\infty f_s(v_\parallel) dv_\parallel=1$,
-
$k_\parallel$ is the wavenumber parallel to the velocity $v_\parallel$ (and also parallel the magnetic field, which otherwise does not feature here; the problem should be considered "1D1V" and $\vec B =0$).
When $f(v_\parallel)=\delta(v_\parallel\pm v_b)$ we need to integrate by parts
$$
1 + \Sigma_s \frac{\Pi_s^2}{\omega}\int_{-\infty}^{\infty} \frac{\partial f_s(v_\parallel)}{\partial v_\parallel}\frac{v_\parallel}{\omega-k_\parallel v_\parallel} dv_\parallel=1 - \Sigma_s \frac{\Pi_s^2}{\omega}\int_{-\infty}^{\infty} f_s(v_\parallel) \frac{\partial }{\partial v_\parallel}\left(\frac{v_\parallel}{\omega-k_\parallel v_\parallel}\right) dv_\parallel
$$
so
$$
\frac{\partial }{\partial v_\parallel}\left(\frac{v_\parallel}{\omega-k_\parallel v_\parallel}\right) = \frac{1}{\omega-k_\parallel v_\parallel}+\frac{k_\parallel v_\parallel}{(\omega-k_\parallel v_\parallel)^2}=\frac{(\omega-k_\parallel v_\parallel+k_\parallel v_\parallel)}{(\omega-k_\parallel v_\parallel)^2}
$$
hence
$$
1 - \Sigma_s \frac{\Pi_s^2}{\omega}\int_{-\infty}^{\infty} f_s(v_\parallel)\frac{\omega}{(\omega-k_\parallel v_\parallel)^2}dv_\parallel=0
$$
Two Dirac delta function species, one with $+v_b$ and another with $-v_b$ gives
$$
1 - \frac{\Pi_s^2}{(\omega-k_\parallel v_b)^2}-\frac{\Pi_s^2}{(\omega+k_\parallel v_b)^2} = 0
$$
In dimensionless units $x=\omega/\Pi_s$ and $u=v_b k_\parallel / \Pi_s$
$$
1 - \frac{1}{(x-u)^2}-\frac{1}{(x+u)^2} = 0
$$
reveals
$$
x=\pm \sqrt{u^2+1\pm\sqrt{4u^2+1}}
$$
For fastest growing mode $\frac{\partial x}{\partial u}=0$ for which we find $u=\frac{\sqrt{3}}{2}$ and there are four solutions for $x$, two real $x=\pm\frac{\sqrt{15}}{2}$ and two imaginary $x=\pm\frac{i}{2}$.
All these solutions share $u=\frac{\sqrt{3}}{2}$ for which $v_b k_{\parallel,max} = \frac{\sqrt{3}}{2}\Pi_s$. If the wavenumber of the fastest growing mode is resolved by the PIC code, the nonlinear self-consistent solution will be dominated by the dynamics of the fastest growing mode and exhibit it's growth rate, $\gamma_{max} = \frac{\Pi_s}{2}$ (where $\gamma$ is typically used in this context to represent the imaginary component of the complex frequency). Note that the unstable solutions have zero real frequency.
Employing periodic boundary conditions and tuning the parameters such that exactly $M$ modes of the fastest growing mode fits in the box in the direction of $v_\parallel$ of length $L$, ( $k_{\parallel,max}=M\frac{2\pi}{L}$), we arrive at:
$$
\gamma_{max} = \frac{\Pi_s}{2} = v_b M\frac{2\pi}{L\sqrt{3}},\\
\Pi_s=v_b M \frac{4\pi}{\sqrt{3}L},\\
n=v_b^2 M^2 \frac{m\epsilon_0}{q^2} \frac{16\pi^2}{3 L^2}
$$
Hence if we initialise all particles as the same species (not two distinct species that differ by sign of drift speed) and randomly assign one of half of them with $+v_b$ and the other half with $-v_b$, we get the case where the particle weight is controlled by the number of particles and the total number density $n_T = 2n$. Hence particle weight $w=\frac{n_T L_x L_y}{N_{comp,part}}$ where $L_x = 1$, $L_y=1$ and $N_{comp,part}$ are the lengths of the domain in $x$ and $y$ (both set to unity) and the total number of computational macro-particles used in the simulation, respectively.
In general the particle weight is
$$
w = \frac{2nL_x L_y}{N_{comp,part}} = v_b^2 M^2 \frac{m\epsilon_0}{q^2} \frac{16\pi^2}{3 L^2}\frac{2L_x L_y}{N_{comp,part}}
$$
to squeeze an integral number of maximally unstable modes into the simulation domain.
So for $v_b = M = m = q = \epsilon_0 = L_x = L_y = 1$
- $\gamma_{max} = \frac{\Pi_s}{2} = \frac{\Pi_T}{2\sqrt{2}}=\frac{2\pi}{\sqrt{3}}$
- $n_T = \frac{32\pi^2}{3}$
- $w=\frac{32\pi^2}{3N_{comp,part}}$
- $\Pi_T = \sqrt{\sum_{s=l,r}\Pi_s^2}=\Pi_s\sqrt{2}$
