Boltzmann equation and equilibrium distribution function

The free-streaming term in stationary 1D Boltzmann equation satisfies the equilibrium distribution function, that is: $$\text{Fr}\{f\} := v(k)\dfrac{\partial f(x,k)}{\partial x} - \dfrac{d\varepsilon_\text{sub}(x)}{dx}\dfrac{\partial f(x,k)}{\partial k}$$ equals zero by putting $f_0(x,k) = \exp\left(-\dfrac{k^2}{2m}-\varepsilon_\text{sub}(x)\right)$, with $v(k)=\dfrac{k}{m}$.

However, after discretization: $$\text{Fr}\{f(x_i,k_j)\} := \dfrac{k_j}{m}\Big[f(x_{i+1},k_j) - f(x_i,k_j)\Big] - \dfrac{\varepsilon_\text{sub}(x_{i+1})-\varepsilon_\text{sub}(x_i)}{\Delta x} \Big[f(x_i,k_{j+1}) - f(x_i,k_j)\Big]$$ the equilibrium function $f_0(x_i,k_j) = \exp\left(-\dfrac{k_j^2}{2m}-\varepsilon_\text{sub}(x_i)\right)$ no longer is a solution to our discretized equation (Of course, the solution gets closer and closer to the equilibrium distribution function by refining the grid).

Since it is a very important property for our discretization to satisfy $\text{Fr}\{f_0\} = 0$ by construction, could you suggest a suitable discretization scheme in this regard?

It might be quite important to know why you need your numerical scheme to produce exact solution for special function $f_0$. I am familiar e.g. with such requirement for so called density driven flows where for special form of given density that is variable only in the direction of gravity one expects that the solution, here a pressure, is approximated (or used exactly) such that no artificial numerical flow is created where the exact flow is zero. Another very well-known case are so called well-balanced numerical schemes for shallow water problems that produce zero numerical flow for hydrostatic situations for arbitrary bottom topography. Therefore it might be useful to understand what is the motivation in your case.
In your case it might be that instead of solving for $f$ you solve for $\tilde f$ such that $f=\tilde f - f_0$ where $f_0$ is your equilibrium distribution that should be uniquely defined for your input data. To obtain the PDE for $\tilde f$ you plug $f$ into your PDE. If your problem is $\hbox{Fr}\{f\}=0$ then you obtain $\hbox{Fr}\{\tilde f\}=\hbox{Fr}\{\tilde f_0\}$, and , clearly, the solution is $\tilde f = f_0$, if you use the same discretization for the both sides (that is natural).
Another solution to your requirement is a special interpolation in the construction of numerical approximation, but this looks quite complicated in your case if the function $\epsilon_{\hbox{sub}}$ can have a general form. Where such procedure is successful is e.g. a numerical solution for one-dimensional stationary advection-diffusion equation where in the case of constant input coefficients one obtains an exponential form of exact solution. Motivated by this, one propose a numerical scheme using exponential profile assumption for numerical solution between two grid points that is appropriate also for variable input coefficients and as a "side product" it can produce the exact solution in the case of constant input coefficients.