Why are $L^2$ functions important? From reading around I have three hypotheses:
All I really know is I keep seeing the term "square-integrable" tossed around as if it's obviously important, obviously the condition one would want to impose, etc. But why?
The reason that spaces of square integrable functions arose in the first place was to study the orthogonal trigonometeric (Fourier) series. Interestingly, Parseval had already noted in 1799 the equality that now bears his name: $$ \frac{1}{\pi}\int_{-\pi}^{\pi}f(x)^{2}\,dx = \frac{1}{2}a_{0}^{2}+\sum_{n=1}^{\infty}a_{n}^{2}+b_{n}^{2}, $$ where $a_{n}$, $b_{n}$ are the (Fourier) coefficients $$ a_{n}=\int_{-\pi}^{\pi}f(x)\cos(nx)\,dx,\;\;\; b_{n}=\int_{-\pi}^{\pi}f(x)\sin(nx)\,dx. $$ This comes out of the orthogonality conditions for the $\sin(nx)$, $\cos(nx)$ terms in the Fourier series. No definite connection was seen between Euclidean N-space and the above at that time; such a connection took decades to evolve. But square-integrable functions gained interest in the early 19th century, and especially after the early 19th century work of Fourier.
It took some time to see a general Cauchy-Schwarz inequality, and to begin to see a connection with geometry, eventually leading to inner-product space abstraction for the space of square-integrable functions. The CS inequality wasn't widely known until after the 1883 publication of Schwarz, even though essentially the same result was published in 1859 by another author. Hilbert proposed his $l^{2}$ space by the early 20th centry as an abstraction of the square-summable Fourier coefficient space, but also a abstraction of finite-dimensional Euclidean space. The connection with square-integrable functions was already firmly established.
In hindsight we can see good reasons that square-integrable functions are connected with energy, and other Physics concepts, but the abstraction seems to have been dictated more out of solving equations using 'orthogonality' conditions. Of course many of the equations arose out of solving physical problems; so it's also hard to separate the two. Now, after the fact, there is interpretation of the integral of the square of a function. On the other hand, the Mathematical abstraction of dealing with functions as points in a space, with distance and geometry on those points has been even more far-reaching, and a great part of the impetus for modern abstract and rigorous Mathematics.
Note: All of this happened before Quantum Mechanics.
Reference: J. Dieudonne, "History of Functional Analysis".
There are several reasons, you named some. You need to work with the space to get a feeling why these things are important.
One reason you did not mention which makes them very popular is the fact that they are very well suited to study elliptic partial differential equations, esp of second order. By partial integration, assuming zero boundary conditions, the Laplacian is intimately connected to the square integral of the gradient
$$ \int \Delta u = \int \langle\nabla u, \nabla u\rangle $$
This makes the integral on the right an important object of study and, naturally, you will look at those functions for which this is finite.
Since the Laplacian is the model operator for elliptic second order differential equations (which may actually be true because it is related to the term in the above equation, which, in a sense, behaves like an energy term) this is one of the very fundamental reasons why these spaces are so important.
In addition, there is a very fundamental isometry of $L^2$, the Fourier transform. This is also intimately connected to the study of elliptic PDE. As an isometry it behaves particularly well on $L^2$ and allows to use the full power of Hilbert space theory to be applied to the study of elliptic PDE.
The definition of $L^p$ functions is motivated by their behavior under the operation of integration: Strictly speaking, the elements of $L^p$ are not functions, but equivalence classes of functions, where two functions are equivalent iff they disagree at most on a set of measure zero---in other words, the equivalence classes keep track only of the functions' behavior under integration and discard finer information.
Now, one can show that the dual space to $L^p$ is canonically isomorphic to $(L^q)^*$, where $\frac{1}{p} + \frac{1}{q} = 1$: Explicitly, we can identify any function $f \in L^p$ with the functional $$g \mapsto \int fg$$ on functions $g \in L^q$.
From this point of view, the space $L^2$ is special in that $p = 2$ is the only value for which $p = q$: In particular, the space is dual to itself with respect to the above bilinear identification $L^2 \cong (L^2)^*$, which by dualizing we can think of as a symmetric bilinear form $L^2 \times L^2 \to \mathbb{R}$. In particular, $\int f^2 > 0$ if $f \neq 0$, and so this bilinear form is actually a (positive definite) inner product. (Notice that the equivalence class formulation is necessary for this inequality, as there are functions nonzero on some nonempty measure zero set, i.e., nonzero functions with zero square integral.) In short, if one wants to work with a class of functions on some space $X$ and (1) is concerned foremost with integrability on $X$ and (2) wants a natural identification of the space with its dual, one settles very quickly in the $L^2$ setting. One gets a good deal of additional structure in this setting for free, including a canonical isometry of $L^2$ with itself, the Fourier transform (this can be formulated in other settings too but it is particularly nice here), but perhaps this is already sufficient advertisement.
This is one explanation (perhaps unsatisfactory) for why this class shows up in QM---in this setting one needs an inner product (self-duality) on the space of wavefunctions in order to formulate probabilities and expected values of actions of operators.
In QM a wave $\psi$ strangely gives the full information of particle and $|\psi|^2$ represents probability for finding that particle (just like $|I|^2$ is probability in EM-theory). When total probability is unity, $||\psi||_2=1$, we have the first connection to $L^2$ functions. Since QM is fundamentally probabilistic, we can only deal with expectations values: $$<f(x)>=\int f(x)|\psi|^2 dx$$ Again, since physical expectations value must be finite, we demand $$\int f(x)|\psi|^2 dx \leqq ||f||_2||\psi||_2=||f||_2 <\infty$$ So again $f\in L^2$ is utmost important for theory.
E.g. for Coulomb potential: $$<V(x)>=<\frac{-ke^2}{|x|}>=-ke^2 \int_{R^3} \frac{1}{|x|}| \psi|^2dx$$
Fourier transform in $L^2$ is important for QM, since that gives us operators. If one considers particle as a wave packet then we know that transform $\hat \psi$ is the probability amplitude for momentum and $|\hat \psi|^2$ is the corresponding probability (since isomorphy $\psi \rightarrow \hat \psi$ ensures that $\hat \psi$ contains the same informations as $\psi$ and due to Parseval's relation $||\hat \psi||_2=||\psi||_2=1$). Now we can calculate those tricky expectations values for momentum $p$: $$<p>=\int_{R^3} p|\hat \psi|^2dp=\int_{R^3} -i\hbar\frac{d}{dx}\psi \centerdot \overline{\psi} dx$$ $$<T>=<\frac{p^2}{2m} >=\int_{R^3} \frac{p^2}{2m}|\hat \psi|^2dp=\int_{R^3} -\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\psi \centerdot \overline{\psi} dx=\frac{{\hbar}^2}{2m}\int_{R^3} |\nabla \psi|^2dx$$
To get some physical results we again need $L^2$ theory:
The lowest expectation value for hydrogen atom, that is the sum of expectation value for kinetic energy and Coulomb potential: $$\text {inf}<T+V(x)>=\text {inf}(\frac{{\hbar}^2}{2m}\int_{R^3} |\nabla \psi|^2dx-ke^2 \int_{R^3} \frac{1}{|x|}| \psi|^2dx)$$
Theorem . If both $f$ and $\nabla f$ belong to space $L^2 (R^3)$ then following inequality holds:$$\int_{R^3} \frac{1}{|x|}| f|^2dx \leqq ||\nabla f||_2 || f||_2$$
[You will have much fun in proving that!]
When you insert that theorem to $\text {inf}<T+V(x)>$ everything boils down to minimizing expression: $$\frac{{\hbar}^2}{2m}||\nabla \psi||^2_2 -ke^2||\nabla \psi||_2$$
Keeping $||\nabla \psi||_2$ as a variable you can use elementary calculus to find the fundamental physical result, the ground state energy for hydrogen: $$\text {inf}<T+V(x)>=-k^2me^4/2\hbar^2.$$ Now you have $||\nabla \psi||_2=kme^2/\hbar^2$ and $\text {inf}<T+V(x)>=-k^2me^4/2\hbar^2$, so you can put them back into expectation value equality: $$-k^2me^4/2\hbar^2=\frac{{\hbar}^2}{2m}\centerdot (kme^2/\hbar^2)^2 - ke^2 \int_{R^3} \frac{1}{|x|}| \psi|^2dx$$ You have $$<\frac{1}{ |x|}>=\int_{R^3} \frac{1}{|x|}| \psi|^2dx=kme^2/\hbar^2=1/a_0$$ That is, the Bohr's radius! $$a_0=\frac{\hbar^2}{kme^2}.$$
Another case from physics is the set of Maxwell's equations in vacuum. There the energy \begin{equation*} \mathcal{E}(t)=\frac{1}{2}\int_{\mathbb{R}^{3}}d\mathbf{x}\{\mathbf{E(x},t% \mathbf{)}^{2}+\mathbf{B(x},t\mathbf{)}^{2}\} \end{equation*} is conserved, i.e. not $t$-dependent, which is easily obtained for smooth field functions $\mathbf{E(x},t\mathbf{)}$ and $\mathbf{B(x},t\mathbf{)}$. Adopting an $L^{2}$-setting this can then be extended to general square integrable fields. Then Maxwell's equations can be given as a unitary evolution on $\mathcal{H}=L^{2}(d\mathbf{x},\mathbb{R}^{3},\mathbb{C}^{6}).$ More general cases such as linear dielectrics can also be dealt with.
