How do I integrate the Poisson equation to determine the electric potential along a particular direction (e.g., $z$)? This question is a sequel of sorts to my earlier (resolved) question about a recent paper.  In the paper, the authors performed molecular dynamics (MD) simulations of parallel-plate supercapacitors, in which liquid resides between the parallel-plate electrodes.  The system has a "slab" geometry, so the authors are only interested in variations of the liquid structure along the $z$ direction.
In my previous question, I asked about how particle number density is computed.  In this question, I would like to ask about how the electric potential is computed, given the charge density distribution.
Recall that in CGS (Gaussian) units, the Poisson equation is 
$$\nabla^2 \Phi = -4\pi \rho$$
where $\Phi$ is the electric potential and $\rho$ is the charge density.  So the charge density $\rho$ is proportional to the Laplacian of the potential.
Now suppose I want to find the potential $\Phi(z)$ along $z$, by integrating the Poisson equation.  How can I do this?
In the paper, on page 254, the authors write down the average charge density $\bar{\rho}_{\alpha}(z)$ at $z$:
$$\bar{\rho}_{\alpha}(z) = A_0^{-1} \int_{-x_0}^{x_0} \int_{-y_0}^{y_0} dx^{\prime} \; dy^{\prime} \; \rho_{\alpha}(x^{\prime}, y^{\prime}, z)$$
where $\rho_{\alpha}(x, y, z)$ is the local charge density arising from the atomic charge distribution of ionic species $\alpha$, $\bar{\rho}_{\alpha}(z)$ is the average charge density at $z$ obtained by averaging $\rho_{\alpha}(x, y, z)$ over $x$ and $y$, and $\sum_{\alpha}$ denotes sum over ionic species.
The authors then integrate the Poisson equation to obtain $\Phi(z)$:
$$\Phi(z) = -4\pi \sum_{\alpha} \int_{-z_0}^z (z - z^{\prime}) \bar{\rho}_{\alpha}(z^{\prime}) \; dz^{\prime} \; \; \; \; \textbf{(eq. 2)}$$
My question is, how do I "integrate the Poisson equation" to obtain equation (2)?  How do I go from  $\nabla^2 \Phi = -4\pi \rho$ to equation (2)?  In paricular, where does the $(z - z^{\prime})$ factor come from? 
Thanks for your time.
 A: I don't know your level of knowledge, so let me start with the very
basic fact that the electric field of a uniformly charged plate is
$$
E=2\pi\sigma,\qquad\left(  1\right)
$$
where $\sigma$ is the surface charge density. To derive this result you can
utilize the Gauss formula:
$$
\Phi=4\pi Q,\qquad\left(  2\right)
$$
where $\Phi$ is the total flux of the electric field through a closed surface
and $Q$ is the total charge in a space bounded by the surface. In the figure
below I depicted charged plate as a blue plane and the closed surface as the
box with green sides.

The flux is only non zero for these green rectangles
$\Phi=2ES$, where $S$ is the area of the rectangles. The total charge inside
the box is $Q=S\sigma$ hence
$$
2ES=4\pi S\sigma\quad\Rightarrow\quad E=2\pi\sigma.
$$
Let's now approximate your system as the set of of plates with surface charge
density $\sigma=\rho\left(  z\right)  \,dz$ where $\rho\left(  z\right)  $ is
the $xy$-averaged charge density. Therefore, the total electric field in a
point $z$ is the difference of the contributions of planes before $z$ and
after $z$ (see figure below):
$$
E\left(  z\right)  =E_{1}\left(  z\right)  -E_{2}\left(  z\right)  ,\qquad(3)
$$
where
$$
E_{1}\left(  z\right)  =2\pi\int_{-z_{0}}^{z}\rho\left(  z^{\prime\prime
}\right)  \,dz^{\prime\prime},\qquad E_{2}\left(  z\right)  =2\pi\int
_{z}^{z_{0}}\rho\left(  z^{\prime\prime}\right)  \,dz^{\prime\prime}.
$$

Thus, the potential $\phi\left(  z\right)  $ has the form:
$$
\phi\left(  z\right)  =-\int_{-z_{0}}^{z}dz^{\prime}E\left(  z^{\prime
}\right)  ,\qquad(4)
$$
with the boundary value $\phi\left(  -z_{0}\right)  =0$. The expression (4) is
the potential required. Let's now simplify it. First of all, I simplify the
expression for the field:
$$
E\left(  z\right)  =2\pi\int_{-z_{0}}^{z}\rho\left(  z^{\prime\prime}\right)
\,dz^{\prime\prime}-2\pi\int_{z}^{z_{0}}\rho\left(  z^{\prime\prime}\right)
\,dz^{\prime\prime}=4\pi\int_{-z_{0}}^{z}\rho\left(  z^{\prime\prime}\right)
\,dz^{\prime\prime}-2\pi\int_{-z_{0}}^{z_{0}}\rho\left(  z^{\prime\prime}\right)
\,dz^{\prime\prime}.
$$
Therefore the potential takes the form:
$$
\phi\left(  z\right)  =-\int_{-z_{0}}^{z}dz^{\prime}E\left(  z^{\prime
}\right)  =-4\pi\int_{-z_{0}}^{z}dz^{\prime}\int_{-z_{0}}^{z^{\prime}}
\rho\left(  z^{\prime\prime}\right)  \,dz^{\prime\prime}-2\pi\left(
z+z_{0}\right)  \int_{-z_{0}}^{z_{0}}\rho\left(  z^{\prime}\right)
\,dz^{\prime}.
$$
To simplify the first term I change the order of integrations (integration
domain is presented in the figure below):
$$
\int_{-z_{0}}^{z}dz^{\prime}\int_{-z_{0}}^{z^{\prime}}\rho\left(
z^{\prime\prime}\right)  \,dz^{\prime\prime}=\int_{-z_{0}}^{z}dz^{\prime
\prime}\int_{z^{\prime\prime}}^{z}\rho\left(  z^{\prime\prime}\right)
\,dz^{\prime}=\int_{-z_{0}}^{z}\left(  z-z^{\prime\prime}\right)  \rho\left(
z^{\prime\prime}\right)  dz^{\prime\prime}.
$$

Finally, we obtain the following result for the potential:
$$
\phi\left(  z\right)  =-4\pi\int_{-z_{0}}^{z}\left(  z-z^{\prime}\right)
\rho\left(  z^{\prime}\right)  dz^{\prime}-2\pi\left(  z+z_{0}\right)
\int_{-z_{0}}^{z_{0}}\rho\left(  z^{\prime}\right)  \,dz^{\prime}.
$$
One can see that the result you presented is valid only for a neutral liquid:
$$
\int_{-z_{0}}^{z_{0}}\rho\left(  z^{\prime}\right)  \,dz^{\prime}=0.
$$
A: I) OP notes that  $\bar{\rho}_{\alpha}$ is the $xy$-averaged charge density
$$
\bar{\rho}_{\alpha}(z) ~:=~ A_0^{-1} \int_{-x_0}^{x_0} \int_{-y_0}^{y_0} dx^{\prime} \; dy^{\prime} \; \rho_{\alpha}(x^{\prime}, y^{\prime}, z), \qquad A_0~:=~4x_0y_0.$$
II) Similarly, eq. (2) must be for the $xy$-averaged potential 
$$ \bar{\Phi}(z) ~:=~ A_0^{-1} \int_{-x_0}^{x_0} \int_{-y_0}^{y_0} dx^{\prime} \; dy^{\prime} \; \Phi(x^{\prime}, y^{\prime}, z),$$
$$\tag{eq. 2}\bar{\Phi}(z) ~=~ -4\pi \sum_{\alpha} \int_{-z_0}^z (z - z^{\prime}) \bar{\rho}_{\alpha}(z^{\prime}) \; dz^{\prime} .$$
Moreover, eq. (2) satisfies the boundary condition 
$$\bar{\Phi}(z=-z_0)~=~0.$$ 
Now differentiate eq. (2) twice wrt. $z$,
$$ \bar{\Phi}^{\prime\prime}(z) ~=~  -4\pi \sum_{\alpha}  \bar{\rho}_{\alpha}(z).$$
This is consistent with the Poisson equation 
$$ \nabla^2 \Phi ~=~ -4\pi \rho, \qquad \rho~=~\sum_{\alpha}\rho_{\alpha},$$
if there is no net electric flux going through the $4$ walls of the box that are parallel to the $z$-directions. 
