Boltzmann distribution of electrons in confining potential I have a particle simulation wherein many non-interacting electrons are trapped in a electric potential well. I would expect, and therefore I initialize according to this, that the electrons would take the shape of a Boltzmann distribution. However, when the simulation is allowed to run for a long time, the charge density converges to a slightly different shape:

In this chart, the horizontal axis denotes position. The red line is the negative of the confining potential. I flipped it just for readability. This is on one vertical axis. On the other vertical axis is both the charge density as the system is initialized (green), and the charge density that the system settles into (blue). The colors are somewhat hard to see, so to clarify, the initialized (green) state is the one that is more negative in the region from around z=20cm to z=60cm
Any ideas as to why the Boltzmann distribution is not the correct initialization?
UPDATE: This is the graph with the electron-electron interaction turned on. It has pretty much settled into equilibrium by now

 A: You can use a Boltzmann distribution at any point in time to describe the density of electrons as a function of position (or potential) only if there are e-e collisions.  These collisions thermalize the distribution of electrons and allow you to define a temperature.  
The distribution you start off with (I assume this is what you mean by "initialization") is an initial condition, and depends on whether electrons had enough time to thermalize before the potential was applied. There is no "right" or "wrong" initialization.
The noise that you're seeing in your simulation looks a lot like a numerical instability or convergence issue.  
Note that you can solve this problem analytically using the electron fluid equation as long as you assume that the a temperature is well-defined at all times (and thus, assume that e-e collisions are always taking place on a fast timescale).
A: Longitudinal averaging
As pointed in the answer by @user3814483 that, in order to obtain a Boltzmann distribution, one needs to assume that there is interaction between particles. While the final formulas of statistical physics are supposedly for non-interacting systems (like the ideal gas), these residual interactions are always assumed to exist, since it is these interactions that lead to the establishment of the thermodynamic equilibrium. In theoretical arguments one may assume these interactions to be very small, which means that we have to wait for a very long time for the equilibrium to be reached. In simulation the strength of the interactions and the duration of the simulation have to be balanced.
ensemble averaging
Another potential approach is using the ensemble average - that is running the simulation multiple number of times, but using different initial conditions (all drawn, e.f., from the Boltzmann distribution). We then expect that the average result of the multiple simulations would resemble the Boltzmann distribution.
Remark: the equivalents of longitudinal average with interactions and the ensemble averaging with randomized initial conditions is what we call ergodicity.
A: To get the Boltzmann distribution, you suppose not only that the potential does not change on the size scale of the contact-interaction, but also that speed of colliding particles are independent. This cannot hold if the material is dense, i.e. if the distance between collisions are in the same order of magnitude as the particle-particle distance.
Is you simulation in 2D or in 3D? How many particles you have? In what volume? How large they are?
