In the quantum harmonic oscillator problem, how would one go about calculating
$$\left\langle n\left|\frac{1}{X^2}\right|n\right\rangle$$
using raising and lowering operators $a^{\dagger}, a$ only, where $X\propto a + a^{\dagger}$ is a linear combination of the raising and lowering operators?
It would also be helpful if someone could refer me to materials on properties of functions of quantum operators.
The operator $X^{-2}$ does exist and is self-adjoint as it follows from standard spectral theory. Its domain is $$D(X^{-2}) := \left\{\psi \in L^2(\mathbb R, dx) \:\left|\: \int_{\mathbb R} x^{-4} |\psi(x)|^2 dx \right.< +\infty\right\}$$ and thereon $$(X^{-2}\psi)(x) := x^{-2}\psi(x)\:.$$ The eigenfunctions of the Hamiltonian operator of the harmonic oscillator are of the form $\psi_n(x) = H_n(x) e^{-x^2/2}$ (with ``normalized'' values of the physical quantities, $m,\omega, \hbar$ appearing in the formula of Hamiltonian), where $H_n$ is a polynomial of degree $n$. Therefore only polynomials $H_n$ which can be factorized as $x^2Q_n(x)$ where $Q_n(x)$ is a polynomial of degree $n-2$ define elements of $D(X^{-2})$. The $H_n$ are the the well-known Hermite polynomials. It is known that $H_{2n}(0) \neq 0$ and $H_{2n-1}(x)$ tends to $0$ with the same order as $x$ for $x\to 0$. Therefore no $\psi_n$ belongs to $D(X^{-2})$ and $\langle \psi_n|X^{-2}\psi_n\rangle$ does not exist.
The use of the operators $a$ and $a^\dagger$ is delicate and, in a sense, dangerous for these type of mathematical problems, because identities like $$X^n= c^n(a+a^\dagger)^n$$ only hold in a subdomain $D$ of $D(X^n)$ given by the (dense) finite span of all functions $\psi_n$, even if the restriction of $X^n$ to that subdomain completely determins $X^n$ itself (as $D$ is a core of $X^n$).
However, I suspect that you are interested in $\langle X^{-2} \rangle_{\psi_n}$,i.e., the expectation value of $X^{-2}$ in the state representaed by $\psi_n$.
To this end it is worth stressing that, if $\psi \in D(A)$ (and I henceforth assume that $\psi$ is normalized) then $$\langle A\rangle_\psi = \langle \psi| A \psi\rangle\:.\tag{1}$$ However, the general definition of $\langle A\rangle_\psi$ does not require $\psi \in D(A)$, but $\psi \in D(\sqrt{|A|})$ is sufficient and in this case $$\langle A \rangle_\psi := \int_{\sigma(A)} \lambda d\mu^{(A)}_{\psi}(\lambda)\tag{2}\:,$$ where $\mu^{(A)}_{\psi}(E) := \langle \psi| P^{(A)}(E)\psi \rangle$ with $E\subset \sigma(A) \subset \mathbb R$ a Borel set and $P^{(A)}$ the spectral measure of $A$, so that $A = \int_{\sigma(A)} \lambda dP^{(A)}(\lambda)$. If $\psi \in D(A)$ which is a subset of $D(\sqrt{|A|})$, it turns out that, as theoretical physicists assume from scratch, (1) holds true, howerver the true definition of $\langle A \rangle_\psi$ is (2).
In the considered case where $A= X^{-2}$, using the given definitions one finds $$\langle X^{-2} \rangle_\psi := \int_{\mathbb R} x^{-2} |\psi(x)|^2 dx$$ provided the integrand is absolutely integrable which means $\psi \in D(\sqrt{|X^{-2}|})$. This is the case for $n= 2k+1$ when $\psi=\psi_n$ because $x^{-2} |\psi_{2k+1}(x)|^2$ is bounded in a neighborood of $x=0$ as said above.
Assuming it, one may go on formally. For instance (with $c$ as above) $$\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}$$ The computation of $\langle X^{-2} \rangle_{\psi_{2k+1}}$ can be arranged similarly starting from $$\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$$ and taking advantage of CCR.
Ali Moh's answer is nice. Here is a different perspective.
Let me just completely forget about creation and annhilation operators and try to construct this matrix element using the wave function for the harmonic oscillator. I'll interpret $X$ as the position operator (since that seems like what you have in mind).
Then let's take a look at this expectation value in the ground state \begin{equation} \mathcal{E}_0 = \langle 0 | \frac{1}{X^2} | 0 \rangle \end{equation} The ground state wave function is $\langle x | 0 \rangle = \psi(x) = \pi^{-1/4} e^{-x^2/2}$. So \begin{equation} \mathcal{E}_0 = \int_{-\infty}^{\infty} dx \psi^*(x) \frac{1}{x^2} \psi(x) = \frac{2}{\sqrt{\pi}} \int_0^\infty dx \frac{e^{-x^2}}{x^2} \end{equation} This integral diverges. To see this we can put in a cutoff on the lower bound: \begin{equation} \mathcal{E}_0[\epsilon] = \frac{2}{\sqrt{\pi}}\int_\epsilon^\infty \frac{e^{-x^2}}{x^2} = 2\left(\frac{1}{\epsilon} - \sqrt{\pi} + O(\epsilon)\right) \end{equation} To me this looks related to the problem that Ali Moh identified, in that $|x\rangle$ is a zero eigenvector of $X$ so we have to be careful near $x=0$. Empirically (ie, trying out the first 10 or so on mathematica), it looks like the integral diverges for all even $n$. (The odd $n$ expectation values are just zero, because the wave function is odd in that case).
An alternative approach is to regulate your operator in such a way that it is inverible. That is, instead of inverting $X^2$, let's invert $X^2 + \ell^2$, where $\ell$ is a constant with dimensions of length.
Then \begin{equation} \mathcal{E}_0 = \frac{2}{\sqrt{\pi}}\int_0^\infty \frac{e^{-x^2}}{x^2 + \ell^2} = \frac{\sqrt{\pi}}{\ell} + \cdots \end{equation} This is just a different way to regularize the integral, but I think in this form it is more clear that the problem is with the operator $1/X^2$ being singular.
when we define $X=a+ a^\dagger$, $\frac{1}{X}$ and any power thereof do not exist.
Proof:
Consider the state $\left|x\right\rangle\propto e^{-a^{\dagger 2} + 2 x\, a^\dagger}\left| 0\right\rangle$, you can show that $X\left|x\right\rangle = x \left|x\right\rangle$, and in particular if you choose $x=0$, then $X\left|x=0\right\rangle = 0$.
Now from elementary linear algebra, an operator with a zero eigenvalue is singular and non-invertible, that is $X^{-1}$ does not exist.
edit: if we gloss over the abstract question of invertibility (because of the subtleties of the infinite hilbert space), you can look at $$X^{-1}X=1$$ \begin{align} \Rightarrow X^{-1}&\left|1\right\rangle = \left|0\right\rangle\\ X^{-1}&\left(\left|0\right\rangle+\sqrt{2}\left|2\right\rangle \right) = \left|1\right\rangle \\ X^{-1}&\left(\sqrt{2}\left|1\right\rangle+\sqrt{3}\left|3\right\rangle \right) = \left|2\right\rangle \\ \ldots \end{align} to convince your self that these conditions are mutually incompatible, in the sense that you cannot write an expansion in powers of $a$ and $a^\dagger$ that would satisfy this conditions. (equivalently write down $\left\langle n |X|n'\right\rangle$ in matrix form and convince yourself you cannot multiply it by any matrix that gives unity)
