El operador $X^{-2}$ existe y es auto-adjunto como se desprende de la norma espectral de la teoría. Su dominio es $$D(X^{-2}) := \left\{\psi \in L^2(\mathbb R, dx) \:\left|\: \int_{\mathbb R} x^{-4} |\psi(x)|^2 dx \right.< +\infty\right\}$$ y sobre la misma
$$(X^{-2}\psi)(x) := x^{-2}\psi(x)\:.$$
Las funciones propias del operador Hamiltoniano del oscilador armónico son de la forma $\psi_n(x) = H_n(x) e^{-x^2/2}$ (con `normalizado" los valores de las magnitudes físicas, $m,\omega, \hbar$ que aparecen en la fórmula de Hamilton), donde $H_n$ es un polinomio de grado $n$. Por lo tanto, sólo polinomios $H_n$ que pueden ser factorizados como $x^2Q_n(x)$ donde $Q_n(x)$ es un polinomio de grado $n-2$ definir los elementos de $D(X^{-2})$. El $H_n$ son bien conocidos de los polinomios de Hermite. Se sabe que $H_{2n}(0) \neq 0$ $H_{2n-1}(x)$ tiende a $0$ con el mismo orden como $x$$x\to 0$. Por lo tanto, no $\psi_n$ pertenece a $D(X^{-2})$ $\langle \psi_n|X^{-2}\psi_n\rangle$ no existe.
El uso de los operadores de $a$ $a^\dagger$ es delicada y, en un sentido, peligroso para este tipo de problemas matemáticos, debido a que las identidades como $$X^n= c^n(a+a^\dagger)^n$$
sólo en un subdominio $D$ $D(X^n)$ dado por la (denso) finita del espacio de todas las funciones $\psi_n$, incluso si la restricción de $X^n$ a dicho subdominio completamente determins $X^n$ sí (como $D$ es un núcleo de $X^n$).
Sin embargo, sospecho que usted está interesado en $\langle X^{-2} \rangle_{\psi_n}$,es decir, la expectativa de valor de $X^{-2}$ en el estado representaed por $\psi_n$.
Para este fin, es importante destacar que, si $\psi \in D(A)$ (y yo de ahora en adelante suponer que $\psi$ está normalizado), a continuación, $$\langle A\rangle_\psi = \langle \psi| A \psi\rangle\:.\tag{1}$$
Sin embargo, la definición general de $\langle A\rangle_\psi$ no requieren $\psi \in D(A)$, pero $\psi \in D(\sqrt{|A|})$ es suficiente y en este caso
$$\langle A \rangle_\psi := \int_{\sigma(A)} \lambda d\mu^{(A)}_{\psi}(\lambda)\tag{2}\:,$$
donde $\mu^{(A)}_{\psi}(E) := \langle \psi| P^{(A)}(E)\psi \rangle$
con $E\subset \sigma(A) \subset \mathbb R$ un conjunto de Borel y $P^{(A)}$ el espectral medida de $A$, por lo que el $A = \int_{\sigma(A)} \lambda dP^{(A)}(\lambda)$. Si $\psi \in D(A)$ que es un subconjunto de a $D(\sqrt{|A|})$, resulta que, como los físicos teóricos que asumen a partir de cero, (1) es verdadera,
howerver la verdadera definición de la $\langle A \rangle_\psi$ (2).
En el caso de que $A= X^{-2}$, utilizando la definiciones que uno encuentra
$$\langle X^{-2} \rangle_\psi := \int_{\mathbb R} x^{-2} |\psi(x)|^2 dx$$
siempre el integrando es absolutamente integrable lo que significa $\psi \in D(\sqrt{|X^{-2}|})$. Este es el caso de $n= 2k+1$ al $\psi=\psi_n$ porque $x^{-2} |\psi_{2k+1}(x)|^2$ está delimitada en un neighborood de $x=0$ como se dijo anteriormente.
Suponiendo que, uno puede ir en formalmente. Por ejemplo (con $c$ anterior)
$$\langle X^{-2} \rangle_{\psi_{1}} = c^{-2} \langle \psi_0| a \frac{1}{(a+a^\dagger)^2} a^\dagger \psi_0 \rangle = c^{-2} \langle \psi_0| (a+ a^\dagger) \frac{1}{(a+ a^\dagger)^2} (a+ a^\dagger) \psi_0 \rangle = c^{-2} \langle \psi_0|\psi_0\rangle =c^{-2}$$
El cálculo de $\langle X^{-2} \rangle_{\psi_{2k+1}}$ puede ser organizado de manera similar a partir de
$$\langle X^{-2} \rangle_{\psi_{2n+1}} = \frac{c^{-2}}{2n+1} \langle \psi_0| (a+a^\dagger) a^{2n} \frac{1}{(a+a^\dagger)^2} (a^{\dagger})^{2n}(a+a^\dagger)\psi_0\rangle$$
y tomar ventaja de CCR.