I am interested in solving a generalized eigenvalue problem to find the natural sloshing frequencies of an arbitrarily shaped basin. The physics is described with a non-linear boundary condition that can be linearised to a simple Robin type boundary condition.

I can find very little information on how to handle Robin type boundary conditions for Galerkin BEM.
My (probably quite naive) idea was to impose the condition in strong form, simply replacing the normal derivative of the potential by the appropriate scalar multiple of the potential itself in the representation formula:

I would really appreciate any suggestions on the correctness of this approach. Best, Duncan

This is a correct approach to handle the Robin boundary condition. You can substitute the boundary condition directly into the representation formula to eliminate either the Dirichlet or Neumann data. Then, you can take the Dirichlet or Neumann trace at both sides to obtain a boundary integral formulation.

A different approach is to use the Dirichlet-to-Neumann or Neumann-to-Dirichlet map directly on the Robin boundary condition. This will lead to the same boundary integral formulation as before. This approach is justified in section 7.4 of Steinbach, 2008 (https://doi.org/10.1007/978-0-387-68805-3).

Notice that with the Robin boundary condition you use, the right-hand side will be zero. Also, I assumed you use an interior formulation, since for the exterior formulation the left-hand side of the representation formula is the scattered field (not the total field).