The justification is a standard proof. You might have come across the potential of a dipole ##\phi(\mathbf{r}) = \dfrac{\mathbf{p} \cdot \mathbf{r}}{4\pi \epsilon_0 r^3}##. Therefore if a material has a dipole moment per unit volume ##\mathbf{P} = \dfrac{d\mathbf{p}}{dV}##, the dipole potential is\begin{align*}
\phi(\mathbf{r}) &= \dfrac{1}{4\pi \epsilon_0} \int_{\mathscr{V}} \dfrac{\mathbf{P}(\mathbf{R}) \cdot (\mathbf{r} - \mathbf{R})}{|\mathbf{r} - \mathbf{R}|^3} d^3 \mathbf{R}
\end{align*}where ##\mathbf{R}## is the integration variable over the source charge distribution (contained with in ##\mathscr{V}##). Letting ##\nabla_{\mathbf{R}} = \left( \frac{\partial}{\partial R_x}, \frac{\partial}{\partial R_y}, \frac{\partial}{\partial R_z} \right)## denote differentiation with respect to the coordinates of ##\mathbf{R}##, notice now that\begin{align*}
\nabla_{\mathbf{R}} \cdot \dfrac{\mathbf{P}(\mathbf{R})}{|\mathbf{r} - \mathbf{R}|} &= \sum_i \dfrac{\partial}{\partial R_i} \dfrac{P_i(\mathbf{R})}{|\mathbf{r} - \mathbf{R}|} \\
&= \sum_i \left[ \dfrac{1}{|\mathbf{r} - \mathbf{R}|} \dfrac{\partial P_i}{\partial R_i} + P_i \dfrac{\partial}{\partial R_i} \dfrac{1}{|\mathbf{r} - \mathbf{R}|} \right] \\
&= \sum_i \left[ \dfrac{1}{|\mathbf{r} - \mathbf{R}|} \dfrac{\partial P_i}{\partial R_i} + P_i \dfrac{(\mathbf{r} - \mathbf{R})_i}{|\mathbf{r} - \mathbf{R}|^3} \right] \\
&= \dfrac{\nabla_{\mathbf{R}} \cdot \mathbf{P}}{|\mathbf{r} - \mathbf{R}|} + \dfrac{\mathbf{P} \cdot (\mathbf{r} - \mathbf{R})}{|\mathbf{r} - \mathbf{R}|^3}
\end{align*}and therefore\begin{align*}
\phi(\mathbf{r}) &= \dfrac{1}{4\pi \epsilon_0} \int_{\mathscr{V}}\left[ \nabla_{\mathbf{R}} \cdot \dfrac{\mathbf{P}}{|\mathbf{r} - \mathbf{R}|} - \dfrac{\nabla_{\mathbf{R}} \cdot \mathbf{P}}{|\mathbf{r} - \mathbf{R}|} \right] d^3 \mathbf{R} \\
&= \dfrac{1}{4\pi \epsilon_0} \int_{\partial \mathscr{V}} \dfrac{\mathbf{P}}{|\mathbf{r} - \mathbf{R}|} \cdot \hat{\mathbf{n}} dS- \dfrac{1}{4\pi \epsilon_0} \int_{\mathscr{V}} \dfrac{\nabla_{\mathbf{R}} \cdot \mathbf{P}}{|\mathbf{r} - \mathbf{R}|} d^3 \mathbf{R} \\
&= \dfrac{1}{4\pi \epsilon_0} \int_{\partial \mathscr{V}} \dfrac{\sigma_{\mathrm{b}}}{|\mathbf{r} - \mathbf{R}|} dS + \dfrac{1}{4\pi \epsilon_0} \int_{\mathscr{V}} \dfrac{\rho_{\mathrm{b}}}{|\mathbf{r} - \mathbf{R}|} d^3 \mathbf{R}
\end{align*}upon defining ##\sigma_{\mathrm{b}} \equiv \mathbf{P} \cdot \hat{\mathbf{n}}## and ##\rho_{\mathrm{b}} \equiv - \nabla_{\mathbf{R}} \cdot \mathbf{P}##.