I am trying to carry out the same relations and I have noticed that the expression (2) is only compatible if the state is acted \begin{align*}a_x^{\dagger}a_y\ket{n_xn_yn_z}\end{align*} with the conjugate complex, that is, by making
\begin{align*} \bra{n_xn_yn_z}a_xa_y^{\dagger}&=\bigg(a_x^{\dagger}a_y\ket{n_xn_yn_z}\bigg)^{\dagger}=\bigg(\sqrt{(n_x+1)n_y}\ket{n_x+1,n_y-1,n_z}\bigg)^{\dagger}\\
=&\sqrt{(n_x+1)n_y}\bra{n_x+1,n_y-1,n_z}
\end{align*}
So
\begin{equation*}
\bra{n_xn_yn_z}a_x^{\dagger}a_y=\bigg(a_xa_y^{\dagger}\ket{n_xn_yn_z}\bigg)^{\dagger}=
\sqrt{n_x(n_y+1)}\bra{n_x-1,n_y+1,n_z}
\end{equation*}
This is the form to find the expression (2).
Later,
\begin{align*}m\braket{100|01m}=-i\braket{010|01m}\hspace{2cm}N=1, n_x=1,n_y=n_z=0\end{align*} and so with the others relations
Starting with the eigenstate $$q=0,l=1,m=1$$ , it holds that the sum of terms is
\begin{equation*}
\ket{011}=\braket{100|011}\ket{100} +\braket{010|011}\ket{010}+\braket{001|011}\ket{001}
\end{equation*}
when $$m=1$$, so the relations are
\begin{align*}
\alpha&=\braket{100|011}=-i\braket{010|011}\\
\beta&=\braket{010|011}=i\braket{100|011}\\
\gamma&=\braket{001|011}=0
\end{align*}
For ortonormality of states , $$\braket{011|011}=1$$, such that
\begin{equation*}
1=|\braket{100|011}|^2+|\braket{010|011}|^2=|\alpha|^2+|\beta|^2
\end{equation*}
Observing,
\begin{align*}
\alpha=-i\beta\\
\beta=i\alpha
\end{align*}
Therefore,
\begin{equation*}
1=|\alpha|^2+|i\alpha|^2=2|\alpha|^2\implies \alpha=\frac{1}{\sqrt{2}}
\end{equation*}
Reemplacing this value of α and β, it is concluded that the state $$\ket{q=0,l=1,m=1}$$ is
\begin{equation*}
\ket{011}=\frac{1}{\sqrt{2}}\ket{100}+\frac
{i}{\sqrt{2}}\ket{010}.
\end{equation*}
I hope it will help. Greetings
