What is the spectrum of $\hat x \hat p + \hat p \hat x$?

This is a nice exercise! It can be completely solved with relatively elementary mathematical techniques.

Let us start by assuming $\hbar:=1$, defining the formally selfadjoint differential operator over smooth functions $$D := \frac{1}{2}(XP +PX) = -i \left(x \frac{d}{dx}+\frac{1}{2}I\right)\:,$$ and proving that it admits a unique selfadjoint extension over natural spaces of functions used in elementary formulation of QM in $L^2(\mathbb{R},dx)$: the space $\cal S(\mathbb{R})$ of Schwartz' functions and ${\cal D}(\mathbb{R}):= C^\infty_c(\mathbb{R})$. Later we will pass to determine the spectrum of the (unique selfadjoint extension of) $D$ we will indicate by $A$.

I stress that without fixing the domain and proving that the operator is selfadjoint thereon or, more weakly, that it admits only one selfadjoint extension on that domain, every physical interpretation as an observable is meaningless and the properties of the spectrum have no clear interpretation.

$\cal S(\mathbb{R})$ and ${\cal D}(\mathbb{R})$ are the most natural and used domains of (essentially) selfadjointness of operators discussed in QM on $L^2(\mathbb{R},dx)$. For instance, the position, momentum, and harmonic Hamiltonian operators are defined thereon giving rise to the known selfadjoint operators.

When $D$ is interpreted as a generator of some one-parameter group of symmetries of a larger group, the domain is fixed in accordance to Garding theory and it may be different form the two cases considered above. Generally speaking, the algebraic properties are not able to fix a selfadjoint extension of the formal observable. For this reason, the analysis of the domain and of selfadjoint extensions is a crucial step of the physical intepretation.

Part 1. To prove that $D$ is essentially selfadjoint, we show that $D$ is the restriction of the selfadjoint generator of a strongly-continuous one-parameter group of unitary operators $U_t$ and we exploit Stone's theorem and a corollary.

If $\psi \in L^2(\mathbb{R},dx)$, we define the natural unitary action of dilation group on wavefunctions $$(U_t\psi)(x):= e^{t/2}\psi(e^tx)\:.\tag{1}$$ (The apparently extra factor $e^{t/2}$ is actually necessary to preserve the norm of wavefunctions.) It is not difficult to prove that, if $t \in \mathbb{R}$, $$\langle U_t \psi|U_t \phi \rangle = \langle \psi|\phi\rangle\:.$$ Furthermore $U_0=I$, $U_tU_s = U_{t+s}$. Finally, it is possible to prove that $$\int_{\mathbb{R}}|e^{t/2}\psi(e^t x) - \psi(x)|^2 dx \to 0 \quad \mbox{for $t\to 0$,}$$ if $\psi \in L^2(\mathbb{R}, dx)$. Stone's theorem implies that there is a selfadjoint operator $A : D(A) \to L^2(\mathbb{R},dx)$ such that $U_t = e^{itA}$ and its dense domain it is just defined by the set of $\psi \in L^2(\mathbb{R}, dx)$ such that, as $t\to 0$, $$\int_{\mathbb R}\left|\frac{(U_t\psi)(x)-\psi(x)}{t}-i \psi'(x) \right|^2 dx \to 0$$ for some $\psi' \in L^2(\mathbb{R},dx)$. Evidently $$A\psi := \psi'$$ Now observe that, if $\psi$ is smooth $$\frac{\partial}{\partial t}|_{t=0} e^{t/2}\psi(e^t x)= i (D\psi)(x)\:.$$ Actually, a suitable use of Lagrange's theorem and Lebesgue's dominate convergence theorem makes stronger the found result to $$\int_{\mathbb R}\left|\frac{e^{t/2}\psi(e^t x)-\psi(x)}{t}-i (D\psi)(x) \right|^2 dx \to 0$$ for either $\psi \in {\cal S}(\mathbb{R})$ or $\psi \in {\cal D}(\mathbb{R})$. Stone's theorem implies that $D$ is a restriction to these dense subspaces of the selfadjoint generator of $U_t$. Actually, since both $U_t {\cal S}(\mathbb{R}) \subset {\cal S}(\mathbb{R})$ and $U_t {\cal D}(\mathbb{R}) \subset {\cal D}(\mathbb{R})$, a known corollary of Stone's theorem implies that $A$ restricted to these spaces admits a unique selfadjoint extension given by $A$ itself.

In other words, $D$ is essentially selfadjoint over ${\cal S}(\mathbb{R})$ and ${\cal D}(\mathbb{R})$ and the unique selfadjoint extensions are exactly the generator $A$ of the unitary group $U_t$ defined in (1).

Part 2. Let us pass to determine the spectrum of $A$. The idea is to reduce to the spectrum of the momentum operator (in two copies) through a (pair of) unitary map(s).

If $\psi \in L^2(\mathbb{R}, dx)$, let us decompose $\psi = \psi_- + \psi_+$, where $\psi_\pm(x) := \psi(x)$ if $x<0$ or $x>0$ respectively, and $\psi_\pm(x) :=0$ in the remaining cases. Evidently $\psi_\pm \in L^2(\mathbb{R}_\pm, dx)$ and the said decomposition realizes the direct orthogonal decomposition $$L^2(\mathbb{R}, dx) = L^2(\mathbb{R}_-, dx) \oplus L^2(\mathbb{R}_+, dx)\:.$$ It is evident form (1) that $U_t L^2(\mathbb{R}_\pm, dx) \subset L^2(\mathbb{R}_\pm, dx)$ so that, also the generator $A$ of $U_t$ admits these orthogonal subspaces as invariant spaces and the spectrum of $A$ is the union of the spectra of the respective restrictions $A_\pm$.

Let us focus on $L^2(\mathbb{R}_\pm, dx)$ defining a unitary map $$V_\pm : L^2(\mathbb{R}_\pm, dx) \ni \psi \mapsto \phi_\pm \in L^2(\mathbb{R}, dy)$$ with $$\phi_\pm(y) = e^{\pm y/2}\psi_+(\pm e^{\pm y})\tag{2}\:.$$ With this definition, it is clear that $$\phi_\pm (y+t) = e^{\pm t/2}e^{\pm y/2}\psi_\pm(\pm e^{\pm t}e^{\pm y})\:,$$ which means $$e^{itP}V_\pm = V_\pm e^{\pm itA_\pm}\:,$$ where $P$ is the standard momentum operator.

Since $V_\pm$ is unitary, $$\sigma_c(A_\pm) = \sigma_c(\pm P) = \mathbb{R}\:,\quad \sigma_p(A_\pm) = \sigma_p(\pm P) = \emptyset\:.$$ We conclude that $$\sigma(A)= \sigma_c(A) = \mathbb{R}\:.$$

The introduced construction also permits us to construct a family of improper eigenvectors of $A$, exploiting the fact that $P$ has a well-known generalized basis of $\delta$-normalized eigenfunctions $$\phi_k(y) = \frac{e^{iky}}{\sqrt{2\pi}}\:, \quad k \in \mathbb{R} \equiv \sigma_c(P)\:.$$ Taking advantage of the unitaries $V_\pm$, we invert (2) and conclude that $A$ admits a generalized basis of $\delta$-normalized eigenfunctions (check computations please, I think this is the first time I do them!) $$\psi^{(k)}_\pm(x) = \frac{(\pm x)^{1\mp ik}}{\sqrt{2\pi}}\quad \mbox{if $x\in \mathbb{R}_\pm$,}\quad \psi^{(k)}_\pm(x) =0 \quad \mbox{otherwise}\:.$$ Notice that, for every $k\in \mathbb{R}$, there is a couple of independent eigenfunctions, so that the spectrum is twice degenerate.

ADDENDUM. The found operator $A$ (unique selfadjoint extension of $D$) is one of the three generators of a unitary representation of the conformal group $PSL(2,\mathbb{R})$, acting on the compactified real line, the one associated to pure dilations. I remember that many years ago I published a paper on the subject, but I do not remember if I analysed the spectrum of $A$ there...

This operator is an element of the $\mathfrak{su}(1,1)$ Lie algebra. This Lie algebra is spanned by $\{\hat K_0,\hat K_\pm\}$ where these are conveniently realized (for your question) in terms of harmonic oscillator operators: $$ \hat K_0=\frac{1}{4}\left(\hat a^\dagger \hat a+\hat a\hat a^\dagger\right) \, ,\quad \hat K_+=\frac{1}{2}a^\dagger a^\dagger\, ,\quad \hat K_-=\frac{1}{2}\hat a\hat a\, . $$ You can easily convert to the combinations of $\hat x$ and $\hat p$.

The spectrum of $\hat x\hat p+\hat p\hat x$ is continuous and its eigenstates are $\delta$-normalized.

See:

Lindblad, G. and Nagel, B., 1970. Continuous bases for unitary irreducible representations of $ SU (1, 1) $. In Annales de l'IHP Physique théorique (Vol. 13, No. 1, pp. 27-56)

and more generally search for SU(1,1) parabolic basis to find literature on this topic.

Note that the construction of eigenstates is quite technical and not necessarily trivial.