1) because the potential energy function has even symmetry, the solutions to Schrodinger's Equation results in purely even or odd solutions. Consider the Time Independent Schrodinger Equation (TISE) below
$$-\frac{\hbar^2}{2m}\nabla^2 \psi(x) + U(x)\psi(x) = E\psi(x)$$
But the potential contains the symmetry that $U(-x) = U(x)$. What this implies is that under the parity transform $$(x \implies -x)$$
$$-\frac{\hbar^2}{2m}\nabla^2 \psi(-x) + U(-x)\psi(-x) = E\psi(-x)$$
Or because of the potential's even symmetry
$$-\frac{\hbar^2}{2m}\nabla^2 \psi(-x) + U(x)\psi(-x) = E\psi(-x)$$
What this means is that both wavefunctions $\Psi(x)$ AND $\Psi(-x)$ solve this particular wave equation. Further, because the Schrodinger Equation is a linear Differential Equation, any linear combination of these two solutions must also be a solution to the equation. therefore, solutions exist such that:
$$\psi_{even}(x) = \psi(x) + \psi(-x)$$
$$\psi_{odd}(x) = \psi(x) - \psi(-x)$$
These follow from the definitions of even and odd functions which is easily verified.
2) Concerning the modulus squared, $\psi^* \psi$, if $\psi$ is always either even or odd, then when multiplying an odd function with itself you get an even function. likewise, an even function multiplied with itself will be even. Consider trivial polynomials for example,
$x^3$ is odd, but $(x^3)(x^3) = (x^3)^2 = x^6$ which is obviously even.
3) Lastly the fact that $\left<x \right> = 0$ follows from the previous statement, that $|\psi|^2$ is always even for potentials with even symmetry, like the SHO. This is because
$$\left<x \right> = \int_{- \infty} ^{\infty} x|\psi|^2dx$$
In the integrand, there is an odd function ($x$) times the definitely even modulus squared. This creates an odd integrand. For any odd integrand with symmetric limits, this integral must be zero.
There is no physical requirement for a wave function to be smooth, or even continuous. At least if we accept the common interpretation that a wavefunction is nothing else than a "probability amplitude". I.e. it represents, when multiplied with its complex conjugate, a probability density. Now a probability density can be discontinuous or even undefined on sets of Lebesgue measure zero. As a matter of fact, the continuous and smooth wavefunction
$$c(x)=e^{-x^2}$$
and the discontinuous wavefunction
\begin{equation}d(x)=\left\{\begin{aligned}&e^{-x^2}&x\neq 0\\&5&x=0 \end{aligned}\right.\end{equation}
give the same probability distribution. In other words, for the purpose of quantum mechanics they are perfectly indistinguishable, and in fact they belong to the same equivalence class of $L^2(\mathbb{R}^d)$ functions. In again different words, the Hilbert space vector $\psi$ on $L^2(\mathbb{R}^d)$ defined by $c(x)$ and $d(x)$ is the same, and therefore it corresponds to the same quantum state.
For the reason above already, continuity is in my opinion meaningless when we deal with wavefunctions. Of course it is always possible (and convenient) to embed, let's say, functions of rapid decrease $\mathscr{S}(\mathbb{R}^d)$ in $L^2(\mathbb{R}^d)$, i.e. considering, by an abuse of notation, $f\in \mathscr{S}(\mathbb{R}^d)\subset L^2(\mathbb{R}^d)$; however we are implicitly considering the equivalence class of almost everywhere equal functions $[f]\in L^2(\mathbb{R}^d)$ to which $f$ belongs, in place of $f$ itself. So let's suppose now that we modifiy OP's assertion to be:
It is unphysical to consider wavefunctions that do not belong to $\mathscr{C}^1_{L^2}(\mathbb{R}^d)=\Bigl\{[f]\in L^2(\mathbb{R}^d), f\in C^1(\mathbb{R}^d)\text{ with }\int_{\mathbb{R}^d}\lvert f(x)\rvert^2 dx<\infty \Bigr\}$.
Now this is a better defined assertion from a mathematical point of view, and it means that we consider physically meaningful only wavefunctions such that there exist a continuous and differentiable representative in its equivalence class. However, also this requirement is not physical. In fact, there are potentials $V(x)$ such that the Hamiltonian $H=-\Delta/2m +V(x)$ is self-adjoint, and $$e^{-iHt}\Bigl[\mathscr{C}^1_{L^2}(\mathbb{R}^d)\Bigr]\nsubseteq \mathscr{C}^1_{L^2}(\mathbb{R}^d)\; .$$
In other words, there are quantum evolutions for which some "smooth" (in the sense just above) wavefunctions become non-smooth as time passes. One relevant concrete example would be the Coulomb potential $V(x)=\pm \frac{1}{\lvert x\rvert}$.
Best Answer
Very interesting question! I'll start by outlining some of the mathematical basics for this question:
You're looking for a bounded state of some potential $V$, i.e. a non-scattering state. Mathemematically, this translates to an $L^2$-integrable eigenfunction of the Schrödinger operator $-\Delta+V$.
By elliptic regularity, for these functions you instantly get what I understand as your condition (3) (The more precise statement would be that $\psi$ lies in the domain of the self-adjoint version of $\hat{p}$). Basically, the argument here is, that $\psi$ has to be two-times differentiable, by rearranging the Schrödinger equation, for your class of potentials you will have $\Delta \psi\in L^2$, by Fourier transform you then get that the first derivative will also be $L^2$. Hence, $\hat{p}$ is well-defined for $\psi$.
The basic idea why most bound states you will come across are exponentially decaying comes from the following idea: Assume that far away from the origin, $V$ is monotone i.e. it doesn't oscillate. This allows us toestimate $V$ from below by a box potential, which implies that a bounded state of $V$ will be dominated by a bounded state of the box potential. Bounded states of box potentials do exponentially decay, hence the state will decay exponentially. This argument can be made explicit using maximum principles for elliptic PDEs, you may up the mathematical details in e.g.
Berezin and Shubin, The Schrödinger equation (Springer 1991).
So from this argumentation, the answer to your question is almost no for potentials which are monotone far outside. By "almost", I mean that there may be such functions at distinguished, but physically irrelevant values of $E$, for example, consider the potential $$ V(x)=\frac{2-6x^2}{(1+x^2)^2} $$which looks like this:
You may now check that $\psi(x)=\frac{2}{\sqrt{\pi}}\frac{1}{1+x^2}$ is a normalized eigenfunction to this potential with eigenvalue 0. The momentum operator is well-defined for this $\psi$ and $\psi$ obviously decays only polynomially for $x\rightarrow\infty$. So, what happened here? If you try to use the "box" argument, you would compare to a box which is completely negative away from the origin (remember, the box estimates the potential from below), so 0 is already a scattering state for the box! However, looking at the potential, you see that this can only be the case for this exact value of $E$ - for even an $\epsilon$ more energy, you will obtain a scattering state since $V\rightarrow 0$ for $x\rightarrow \infty$; and for an $\epsilon$ less, you will get a bound state you can again estimate by a box, hence it decays exponentially. Since you can't prepare a state with an exact energy, this is not physically relevant. Generally speaking, this phenomenon should only occur at $E=\limsup_{|x|\rightarrow\infty}V(x)$, since this will correspond to the lowest possible energy for scattering states.
So, what can happen if we drop the "monotonicity-far-outside"-condition? I think that in this case, it should be possible to obtain the kind of states you are searching for. My attempt on the construction goes as follows: Let $V$ be a collection of box potentials where the boxes have constant height and grow thinner as $x$ gets bigger, e.g. something that looks vaguely like this:
If the infinitely many discontinuities bother you, the behavior I'll describe should be exactly the same for a smoothed version of this potential.
Now, a bound state of this potential would oscillate around $0$ where $V=-0.5$ and decay where $V=+0.5$. The (exponential) decay rate where $V=0.5$ is always the same, by controlling the width of the boxes, you can exactly control how fast your bound state decays, e.g. you can achieve, that each time you pass the positive part of the box, your amplitude goes down at the rate $1/x^2$. The details are probably very technical and fishy, but I think in principle this should work.