  We derive the kinetic theory of fluctuations in physically and numerically
stable particle-in-cell (PIC) simulations of electrostatic plasmas. The
starting point is the single-time correlations at the simulation start between
the statistical fluctuations of weighted densities of macroparticle centers in
the plasma particle phase-space. The single-time correlations at all time steps
and in each spatial grid cell are then determined from the Laplace-Fourier
transforms of the discretized Klimontovich-like equation for the macroparticles
and Maxwell's equations for the fields as computed by modern PIC codes. We
recover the expressions for the electrostatic field and the plasma particle
density fluctuation autocorrelations spectra as well as the kinetic equations
describing the average evolution of PIC-simulated plasma particles, first
derived by Langdon in 1970, using a macroparticle test approach perturbing a
discretized Vlasovian plasma and then averaging the obtained physical quantity
over the initial macroparticle velocity distribution. We generalize and extend
these results to the modern algorithms in PIC codes and using arbitrary
macroparticle weights. Analytical estimates of statistical fluctuations
single-time correlation amplitudes are derived as a function of the plasma
simulation parameters, using the central limit theorem in the limit of a large
number of macroparticles per cell. The theory is then used to analyze the
ensemble averaging technique of PIC simulations where statistical averages are
performed over ensembles of PIC simulations, modeling the same plasma physics
problem but using different statistical realizations of the initial
distribution functions of the macroparticles.
