Why are wave functions required to vanish at infinity?

In fact, it is nothing but a vague and, strictly speaking, false requirement that can be found in some physically minded books (also of very good level).

There are however physical situations where regularity of the used functions implies that they must vanish at infinity. If one is solving the stationary Schroedinger equation and the potential is sufficiently regular, the eigenvectors must be accordingly regular because of known theorems (especially due to Weyl) on elliptic regularity. In some cases regularity plus the requirement that the function belongs to $L^2$ and some control of the asymptotic value of some derivatives (arising from a nice asymptotic shape of the potential) imply that the wavefunction must vanish at infinity.

On the other hand, from a physical viewpoint, it is impossible to prepare a state in space working arbitrarily far from a given position where physical instruments are localized. It is therefore reasonable to assume that physically realisable states are described by wavefunctions vanishing outside a sufficiently large spatial region. This, in turn, implies a corresponding requirement on physically accessibile observables and their realistic realisation (it does not make sense to assume to fill in the universe with detectors). However these are physical requirements that should not be confused with mathematical restrictions.

From a pure mathematical viewpoint, instead, no requirement about a rapid decay at infinity applies on vawefunctions in $L^2$ (it is easy to construct smooth $L^2(\mathbb R)$ functions with larger and larger oscillations as $x\to \infty$ and corresponding strongly unphysical Hamiltonian operators whose these monsters are eigenfunctions). Nor any strong regularity condition applies. Indeed, these functions are defined up to zero-measure sets and all operators representing observables, in order to be properly self adjoint (not simply Hermitian) are the closure of differential operators whose final domains are therefore Sobolev-like spaces: Derivatives are required to exist in a weak sense at most.

The only exceptions are probably eigenfunctions of Hamiltonian operators with sufficiently regular potentials, where both Sobolev's regularity and elliptic regularity results can be applied and outside singularities of the potential, eigenfunctions are $C^2$ (or even smooth) in a proper sense.

Restriction to the Schwartz space may have sense because most operators have self-adjointness domains including that space (which sometimes is also a core of the operators) and also because Schwartz space is dense in $L^2$.

However it turns out to be a too strong restriction also in some elementary cases. Think of an 1D Hamiltonian operator with a potential with some mild discontinuity. It does not admit eigenfunctions of Schwartz type.

Example.

I construct here a true monster without physical meaning, but permitted by mathematical requirements of QM.

Consider a function $\psi \in L^2(\mathbb R)$ costructed this way. It attains the value $0$ everywhere except for every interval $$I_n= \left(n -\frac{1}{2(n^4+1)^{2}},\:\: n+ \frac{1}{2(n^4+1)^{2}}\right) $$ thus of length $(n^4+1)^{-2}$ and centered on every $n\in \mathbb Z$. In these intervals $$\psi(x) = n^2 \quad x\in I_n\:.$$ It is clear that $\int_{\mathbb R} |\psi(x)|^2 dx$ for some constant $C>0$ satisfies $$\int_{\mathbb R} |\psi(x)|^2 dx \leq C \sum_{n=0}^{+\infty} \frac{n^4}{(n^4+1)^2} <+\infty\:.$$ It is clear that this function oscillates with larger and larger oscillation as $|x|\to +\infty$, but these oscillations are sharper and sharper.

Next, smoothly modify $\psi$ in every interval $I_n$ to produce a smooth non-negative function $0\leq \phi(x) \leq \psi(x)$ with $\phi(n)=\psi(n)$ thus preserving finiteness of the integral of $|\phi|^2$.

Finally define $$\Psi(x) := \phi(x) + \frac{1}{1+x^2}\:.$$ It easy to see that, since $\phi(x)\geq 0$, it must be $$\Psi(x) >0\tag{1}$$ and $$\Psi \in L^2(\mathbb R)\cap C^\infty(\mathbb R)$$ because $\Psi$ is sum of two $L^2$ (smooth) functions and $L^2(\mathbb R)$ is a vector space. A possible qualitative shape of $\Psi$ is represented here

Observe that $\Psi$ is an eigenfunction (a proper eigenvector since it belongs to $L^2$), of the Hamiltonian $$H = -\frac{\hbar^2}{2m}\frac{d^2}{dx^2} + U(x)\:, $$ where $$U(x) := \frac{\hbar^2}{2m}\frac{\Psi''(x)}{\Psi(x)}$$ is smooth by construction, in particular, (1) holding.

It holds $$H\Psi = E \Psi$$ with $E=0$. Obviously all this construction is physically meaningless, but it is permitted by general requirements of QM formalism.

This might be able to answer your question in part. Prof. Moretti has a better answer however.

Consider one dimensional problems, on the Hilbert space $L^2(\mathbb{R})$.

Wave functions, as people call them, are usually eigenfunctions of Hamiltonian, i.e. the energy operator. Even otherwise, they are associated with some physical operator $A$ made from $X$ or $P$. Whatever space $A$ is defined on, these wave functions lie in that space.

Unless you define these operators properly, they are not self-adjoint to begin with; forget about their sums or compositions.

From Stone's theorem, we define $X$ and $P$ as generators of appropriate unitary transformations and the theorem spits out a special domain for self adjointness. The domains (and their images under the map) being unidentical, we can't define compositions like $X\circ X$ and $P\circ P$ over them for example.

We however get away by using Schwartz space $S(\mathbb{R})$ as domain for all of them. It is a dense subset of $L^2$ which maps onto itself under the usual action of $X$ and $P$. Then we can define all sorts of sums and compositions on the Schwartz space. However, their self adjointness is not guaranteed even then. But for the best example we always give, the Harmonic oscillator,

$$H = X\circ X + P\circ P$$

is essentially self adjoint. Thus, its unique self adjoint extension $\bar{H} $ has all its eigenfunctions.

Wave functions need not really be eigenfunctions however. We can prepare states like a Gaussian wave packet. I'm not too sure about what can be allowed for them in full rigor.

So to answer your question to best of my ability, indeed, the wave functions, often do lie in the Schwartz space. But people may be able to find some crazy examples otherwise. I am not too familiar with that.