The normal trick is to set it up as an eigenvalue problem, namely to look instead at
$$\Delta u = 0, ~~\frac{\partial u}{\partial n} = \lambda b(x) u(x)~~ \forall x \in \partial D. (1)$$ 
You first establish that this corresponds to an eigenvalue problem  on $L^2(\partial \Omega)$ for a compact operator, thanks to the compact embedding of $H^{1/2}(\partial D)$ into $L^2((\partial D)$. 

Then the usual machinery rolls out. There is a discrete spectrum, and there are eigenvalues. And if it so happens that $\lambda=1$ is one of those, then your problem has non trivial solutions. 

The advantage of this point of view is that the problem is well posed, there is a clear functional analytic setup, and you don't need much on $b$, some integrability (if it is bounded as you liked to consider, it is fine).

It is done in the following way. Given $f \in H^{-1/2}(\partial D)$ such that $\int_{\partial\Omega} f d\sigma= 0$ consider the so called Neumann-to-Dirichlet map
\begin{eqnarray*}
NtD : H^{-1/2}(\partial D)/\mathbb R &\to& H^{-1/2}(\partial D)/\mathbb R \\
f &\to& u|_{\partial D} : \begin{cases} -\Delta u =0& \text{ in }D \\ \partial_n u =f &\text{ on }\partial D \end{cases}
\end{eqnarray*}
the quotient over $\mathbb R$ means that we impose $\int_{\partial\Omega} f d\sigma= \int_{\partial\Omega} u d\sigma= 0.$

The operator $NtD$ is compact : we are going to use this operator on a restricted domain, where the compactness is very easy to see.
 
Restrict this operator to $L^2(\partial D)/ \mathbb R$, and you can write it explicitely in terms of Fourier coefficients. 
Write
$$
f=\sum_{n=1}^\infty a_n \cos n \theta + b_n \sin n \theta, \sum |a_n|^2 + |b_n|^2 <\infty
$$
then 
$$
NtD(f)= \sum_{n=1}^\infty \frac{a_n}{n} \cos n \theta + \frac{b_n}{n} \sin n \theta.
$$
It is easy to see that it is a compact operator from $L^{2}(\partial D)/\mathbb R$ to $L^{2}(\partial D)/\mathbb R$  (the $n^{-1}$ helps) .

Supppose $\| b \|_{L^\infty(\partial D)} <\infty$. Then
\begin{eqnarray*}
B : L^{2}(\partial D)/\mathbb R &\to& L^{2}(\partial D)/\mathbb R \\
f &\to& b f - \frac{1}{\left|\partial D\right|} \int_{\partial D}bf \end{eqnarray*}
is a continuous, bounded operator from $ L^{2}(\partial D)/\mathbb R$ into itself.

So the map
\begin{eqnarray*}
T : L^{2}(\partial D)/\mathbb R &\to& L^{2}(\partial D)/\mathbb R \\
f &\to& B\left(NtD \left(f\right)\right) 
\end{eqnarray*}
is a compact, linear operator. Its spectrum is discrete, it has a maximal eigenvalue etc. and (1) is 
$$
f=\lambda T f.
$$

-----

Regarding regularity of the solutions. Naturally, they are analytic inside $D(0,\rho)$ for any $\rho<1$, this is just interior regularity of harmonic functions. At the boundary, the regularity is dictated by $b$. Indeed, suppose $b$ is piecewise constant. We have $Du= \partial_r u e_r + \frac{1}{r} \partial_\theta u e_\theta= \lambda b u e_r + \frac{1}{r} \partial_\theta u e_\theta $. If $u$ is regular, this formula can be "extended" a little inside. Then you see immediately that you cannot differentiate $u$ a second time: all terms are differentiable, except one, $b$.