The question deals with a particular operator characterized by the commutation relation $[a, a^\dagger] = 1$. The polynomial to be expanded is also of a special form, $(a+a^\dagger)^n$. Also, the whole polynomial is supposed to be applied to the unique state vector $\vert 0\rangle$ that has the property $a \vert 0\rangle = 0$ (which is the null vector).
So this is a rather specialized function. The only thing that needs to be general here is that it should work for any power $n$. The rules in the question are derived from these very special circumstances. What I'm doing below is an approach that is in principle much more general because it would work with small modifications also for other polynomials and other state vectors to which the operators are applied.
You can solve this problem without using replacement rules. I'll just base the calculation purely on linear algebra, knowing what the matrix of the operator $a$ in the basis spanned by $\vert 0\rangle$, $a^{\dagger}\vert 0\rangle$, etc. looks like:
reduceHO[n_] :=
Module[{a = SparseArray[{Band[{1, 2}] -> Sqrt[Range[n]]}, {n, n} + 1]},
Total[(SuperDagger["a"])^Range[0,n]
MatrixPower[a + Transpose[a], n].UnitVector[n + 1, 1]/Sqrt[Range[0, n]!]]]
The argument of reduceHO is the power $n$ in the polynomial $(a+a^\dagger)^n$ to be expanded. Here is a test:
reduceHO[2]
$$\left(\text{a}^{\dagger
}\right)^2+1$$
reduceHO[5]
$$\left(\text{a}^{\dagger
}\right)^5+10
\left(\text{a}^{\dagger
}\right)^3+15 \text{a}^{\dagger
}$$
This is the expanded operator polynomial which is equivalent to the original one, when applied to the state $\vert 0\rangle$. It is just an attempt at a pretty output format - the detailed formatting could be customized further.
This approach is based on converting the operators to matrices, and then using MatrixPower to do the polynomial. The dimension of the matrix just has to be equal to the largest power of the raising operator $a^{\dagger}$ that can appear in the polynomial.
The basis states have normalization factors in them, which give rise to the Sqrt factors in the calculation of the matrix elements. These are divided out again at the end, leaving only the amplitude factors with which the different powers of the (unnormalized) vectors $(a^\dagger)^\nu\vert 0\rangle$ appear. The MatrixPower first yields a vector when applied to the UnitVector representing the ground state $\vert 0\rangle$, which is converted back to a symbolic expression inside the Total which sums up the list created from the amplitude factors and the symbolic powers $(a^\dagger)^\nu$, written as strings so that they remain inert. The last step could be modified depending on what you intend to do with the result.
Explanations
The matrix elements on which the power is based are taken from here. To show what the matrix has to do with the rules of the original question, here is its form for the case $n=6$:
With[{n = 6},
a = Normal@
SparseArray[{Band[{1, 2}] -> Sqrt[Range[n]]}, {n, n} + 1]];
MatrixForm[a]
$$\left(
\begin{array}{ccccccc}
0 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & \sqrt{2} & 0 & 0 & 0 & 0
\\
0 & 0 & 0 & \sqrt{3} & 0 & 0 & 0
\\
0 & 0 & 0 & 0 & 2 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & \sqrt{5} & 0
\\
0 & 0 & 0 & 0 & 0 & 0 & \sqrt{6}
\\
0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{array}
\right)$$
The state $\vert 0\rangle$ is now represented by the canonical unit vector $\{1, 0, 0, 0, 0, 0, 0\}$. This makes it clear that whenever $a$ is directly applied to $\vert 0\rangle$, the result is the null vector $\{0, 0, 0, 0, 0, 0, 0\}$, because the off-diagonal form of the above matrix shifts the entries of a vector one step to the left. It also multiplies the result by scale factors, but they are needed in order to insure the commutation relation:
MatrixForm[a.Transpose[a] - Transpose[a].a]
$$\left(
\begin{array}{ccccccc}
1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & -6 \\
\end{array}
\right)$$
The last entry is not 1 because of the truncation to a finite dimension; the actual commutation relation can only hold in an infinite-dimensional space (in any finite dimension, taking the trace of $[a, a^\dagger] = 1$ on both sides would cause a contradiction). As long as one avoids such truncation effects by choosing the dimension of the matrix large enough, all the desired algebraic properties of $a$ follow from its matrix form above.
The other matrix appearing in the problem is $a^\dagger = a^T$, the Transpose of $a$. Its effect is to shift a vector in the opposite direction (multiplying it by scale factors again). This is why it's possible to predict that a product with more factors $a$ than $a^\dagger$ will map the state $\vert 0\rangle$ onto the null vector: the net shift is to the left.
