The Green function $G_{pt}$ satisfies $$-\nabla(\epsilon\nabla G_{pt})+\kappa^2 G_{pt}=\delta(t-p),$$ where $\nabla$ describes derivatives with respect to $t$. Physically, that is, $G_{pt}$ is the field at the point $t$ due to a point source at the point $p$, under some (unspecified) set of boundary conditions. Multiply the left side of the linearized PDE, as a function of $t$, by $G_{pt}$, and integrate over $\Omega$. Moving the derivatives over to act on the Green function via two integrations by parts gives the surface terms shown, and a remaining volume integral of $\delta(t-p)\, \phi^{\mathrm{int}}(t)$.
This is standard stuff. Sorry if it's obvious.
The game in the paper seems to be that, where the potential is large and quickly varying, i.e., near the cluster of sources that make up the molecule, the sine term in the original, non-linear PDE is negligible, and one has effectively a Poisson equation. (Or maybe the point is that the fluid has different properties on either side of a physical boundary. Guess who's not a chemist.) Farther away, where the potential is small, one can linearize and get a (modified) Helmholtz equation. Imposing $\phi_t^{\mathrm{int}}=\phi_t^{\mathrm{ext}}$ ties the solutions of those two different PDEs together at the boundary separating "near" from "far."
Eqs. (5) and (6) arise by taking the limit as $p$ approaches the boundary from the inside and outside, respectively. (The principal value gets around the singularity in $G_{pt}$, say, at $p = t$ on the surface.) The boundary values of the potential and electric field need to be related by the two equations in the text just before eq. (5). So, one has two integral equations giving the boundary potential at one point $p$ on the surface, now called $f_p$, in terms of its values $f_t$ everywhere on the surface and the values $h_t$ of the (exterior) normal derivative. Solving these coupled integral equations, which one can hope to do numerically by formulating them as a discretized linear algebra problem, gives a consistent set of boundary data. These boundary data then give the solutions everywhere via eqs. (3) and (4).
So, it shouldn't be obvious that $\phi_p^{\mathrm{int}} = \phi_p^{\mathrm{ext}}$ on the surface. Demanding that equality, and the other relation between the normal derivatives, * is equivalent to solving the PDE*, and must be done numerically.