Okay, I checked and you got it mostly right, except the factor $\mu(u)$ in front of the power of $u$ at the very end, in the main equation.
Let's simplify slightly your notations and write $R(\Omega, u)$ and $R(\Omega, u + \varphi)$ for the Rayleigh quotients. I do not compute using a real parameter ; rather, I perform some generalized Taylor exapnsions to get the full derivative of the quotient, which means that $\varphi$ is to be thought of as small itself. Denote by $\|u\|_q^2$ the denominator of $R$.
The numerator of $R(\Omega, u + \varphi)$ expands as
\begin{equation}
\int_{\Omega} |\nabla u |^2 + \alpha \int_{\partial \Omega} |u|^2 + 2 \int_{\Omega} \nabla u \cdot \nabla \varphi + 2 \alpha \int_{\partial \Omega} u \varphi + \mathcal{O}(\varphi^2).
\end{equation}
Its denominator writes
\begin{align} & \|u\|_q^2 \left( 1 + \frac{q}{\|u\|_q^q} \int_{\Omega} |u|^{q-2} u \varphi + \mathcal{O}(\varphi^2) \right)^{\frac{2}{q}} \\ & = \|u\|_q^2 \left( 1 + \frac{2}{\|u\|_q^q} \int_{\Omega} |u|^{q-2} u \varphi + \mathcal{O}(\varphi^2) \right) . \end{align}
Thanks to the usual $(1+x)^{-1} = 1-x + \mathcal{O}(x^2)$, we get :
\begin{equation} R(\Omega, u + \varphi) = R(\Omega, u) - \frac{2}{\|u\|_q^q} R(\Omega,u) \int_{\Omega} |u|^{q-2} u \varphi + \frac{2}{\|u\|_q^2} \left( \int_{\Omega} \nabla u \cdot \nabla \varphi + \alpha \int_{\Omega} u \varphi \right) + \mathcal{O}(\varphi^2). \end{equation}
Using integration by parts, this rewrites
\begin{equation} R(\Omega, u + \varphi) = \lambda - \frac{2}{\|u\|_q^q} \lambda \int_{\Omega} |u|^{q-2}u \varphi + \frac{2}{\|u\|_q^2} \left( - \int_{\Omega} \Delta u \varphi + \int_{\partial \Omega} (\partial_{\nu} u + \alpha u) \varphi\right). \end{equation}
Finally we are done. The boundary condition that you got was the right one, whereas the main equation is now
\begin{equation} - \frac{\Delta u}{\|u\|_q^2} = \lambda \frac{|u|^{q-2}u}{\|u\|_q^q}, \end{equation}
i.e.
\begin{equation} - {\Delta u} = \lambda \frac{|u|^{q-2}u}{\|u\|_q^{q-2}}. \end{equation} Notice that if we set $q=2$ and $\alpha = 0$, we indeed get back to the usual eigenvalue problem for the Neumann laplacian, which is your first case.
I will attempt to elaborate a bit more on why $q^*$ is important. Recall that the Sobolev embedding theorem tells you in particular that $H^1$ embeds into $L^{q^*}$, where $1 = d(\frac 12 - \frac {1}{q^*})$ and $d$ is the dimension of the ambient space. This embedding means in particular that $L^{q^*}$ scales like $H^1$, in that, for a non zero $f \in H^1$ and $n \in \mathbb{N}$, $\frac{\|f(n \cdot)\|_{H^1}}{\|f(n \cdot)\|_{L^{q^*}}}$ does not depend on $n$, whereas for a general exponent $q$ instead of $q^*$, it will. This means we have three cases to look at :
Case 1 : $q > q^*$.
This is the easiest one and it relies solely on the scaling. Up to a translation, we can assume that $0$ is in the interior of $\Omega$. Let $\rho$ be your favorite function in $\mathcal{C}^{\infty}_c(\Omega)$ localised around $0$ and assume for simplicity that its $L^q$ norm is equal to $1$. Define now $\rho_n$ by $\rho_n(x) := n^{\frac dq} \rho (nx)$. (This defines a family of ultra-thin, ultra-high functions around $0$.)
For any $n \in \mathbb{N}$, we have :
- $\|\rho_n \|_{L^q(\Omega)} = 1$ (thanks to the $n^{\frac dq}$ factor),
- $\|\rho_n \|_{L^2(\partial \Omega)} = 0$ (compact support condition),
- $\|\nabla \rho_n\|_{L^2(\Omega)} = \|\nabla \rho\|_{L^2(\Omega)} n^{d(\frac 1q - \frac{1}{q^*})}$.
Gathering all the equalities above, one obtains hat $R(\Omega, \rho_n) \to 0$ as $n \to \infty$. But the only function cancelling the numerator of $R$ is the zero function, which is not admissible. Thus, in this case, $\lambda = 0$ and there is no minimizer.
Case 2 : $2 \leq q < q^*$.
Here, opposite to what happened in Case 1, you will get both a nontrivial $\lambda$ and the existence of a minimizer.
Let $(u_n)$ be a sequence of functions with the following properties :
- $\|u_n\|_{L^q} \equiv 1$ for all $n$,
- $R(\Omega, u_n) \to \lambda$ as $n \to \infty$.
Because $q \geq 2$ and $\Omega$ has finite measure, the first bullet implies that $\|u_n\|_{L^2} \lesssim 1$. From the second bullet, one has in particular that $(\nabla u_n)$ is bounded in $L^2$. Thus, there exists some $v \in H^1(\Omega)$ such that, up to extraction, $u_n \to v$ strongly in $L^2(\Omega)$ and weakly in $H^1(\Omega)$. By interpolation between $L^2$ and $L^{q^*}$, one also has $u_n \to v$ strongly in $L^q(\Omega)$. In particular, $\|v\|_{L^q(\Omega)} = 1$.
On the other hand, by interpolation between $L^2$ and $H^1$, $u_n \to v$ strongly in $\dot{H}^{\frac34}(\Omega)$ (say). Because the trace operator is continous from $\dot{H}^{\frac 34}(\Omega)$ to $\dot{H}^{\frac 14}(\partial \Omega)$, it follows that $u_n \to v$ strongly in $\dot{H}^{\frac14}(\partial \Omega)$. Also, because $(u_n)$ is bounded in $L^2(\partial \Omega)$, up to extraction, $u_n$ converges weakly in $L^2(\partial \Omega)$ to some function $g$. Identifying distributional limits, one obtains that $g$ equals (the trace on $\partial \Omega$ of) $v$.
Eventually, using Fatou lemma, one has
\begin{equation} R(\Omega, v) \leq \liminf_{n \to + \infty} R(\Omega, u_n) = \lambda. \end{equation}
Thus, $v$ is a minimizer. In particular, since $v$ cannot be $0$ (having for instance unit $L^q$ norm), $\lambda$ is nonzero.
The third case (namely $q = q^*$) seems more involved, I'll try later on. If you want any clarification, please go ahead.