Solution (source code)

= Solution

Set $r_0=q\mathcal E/(m\omega^2)$ and $y=r-r_0$. Completing the square in the <potential energy> gives
$$
\frac12m\omega^2|r|^2-q\mathcal E\cdot r
=\frac12m\omega^2|y|^2-\frac{q^2|\mathcal E|^2}{2m\omega^2}.
$$
Translation leaves the <Laplacian> unchanged, and the shifted three-dimensional <quantum harmonic oscillator> separates into three one-dimensional oscillators. Their energies add, giving
$$
\boxed{E_N=(N+\tfrac32)\hbar\omega-\frac{q^2|\mathcal E|^2}{2m\omega^2},\qquad N=n_1+n_2+n_3.}
$$
The <energy-level degeneracy> is the dimension of the energy eigenspace, or number of linearly independent states with that energy. Counting nonnegative integer triples of total $N$ by <stars and bars> gives
$$
\boxed{d_N=\binom{N+2}{2}=\frac{(N+1)(N+2)}2.}
$$
The constant electric <field> displaces the oscillator and shifts all levels by the same amount, preserving their degeneracies.