Why are wavefunctions with same principal quantum number and different azimuthal quantum number orthogonal
The mathematical answer to this specific question lies in the fact that s and p orbitals are simultaneous, non-degenerate eigenfunctions of the quantum mechanical $\hat{L}^2$ operator, which is hermitian:

$$\begin{align}
\hat{L}^2\left|\mathrm{s}\right&gt; = [l(l+1)\hbar]^2 \left|\mathrm{s}\right&gt; &amp;= 0 \\
\hat{L}^2\left|\mathrm{p}\right&gt; = [l(l+1)\hbar]^2 \left|\mathrm{p}\right&gt; &amp;= 2\hbar^2 \left|\mathrm{p}\right&gt;
\end{align}$$

since $l = 0$ for s orbitals and $l = 1$ for p orbitals.

It is a simple proof that the eigenstates of hermitian operators that correspond to different eigenvalues are necessarily orthogonal. Suppose, for a hermitian operator $A$, we have $A\left|m\right&gt; = a_m\left|m\right&gt;$ and $A\left|n\right&gt; = a_n\left|n\right&gt;$ with $a_m \neq a_n$. Then:

$$\begin{align}
\left\langle m \middle| A \middle| n\right\rangle &amp;= \left\langle n \middle| A \middle| m\right\rangle^* &amp; &amp;\text{(definition of hermiticity)} \\
a_n\left\langle m \middle| n\right\rangle &amp;= (a_m\left\langle n \middle| m\right\rangle)^* &amp; &amp;\text{(extracting the eigenvalues)} \\
&amp;= a_m\left\langle n \middle| m\right\rangle^* &amp; &amp;\text{(eigenvalues of hermitian operator are real)} \\
&amp;= a_m\left\langle m \middle| n\right\rangle &amp; &amp;\text{(definition of inner product)} \\
(a_n - a_m)\left\langle m \middle| n\right\rangle &amp;= 0 \\
\left\langle m \middle| n\right\rangle &amp;= 0 &amp; &amp;(a_m \neq a_n)
\end{align}$$

You can likewise show that any pair of atomic orbitals in the hydrogen atom are orthogonal, using whichever of the three operators $\hat{H}$, $\hat{L}^2$, $\hat{L}_z$ that you need. All atomic orbitals are simultaneous eigenfunctions of these three operators, and have different eigenvalues which correspond to the quantum numbers $n$, $l$, and $m_l$:

$$\begin{align}
\hat{H}\left|n,l,m_l\right&gt; &amp;= E_n\left|n,l,m_l\right&gt; = -\frac{m_e e^4}{2(4\pi\varepsilon_0)^2 \hbar^2 n^2}\left|n,l,m_l\right&gt; \\
\hat{L}^2\left|n,l,m_l\right&gt; &amp;= L^2\left|n,l,m_l\right&gt; = [l(l+1)\hbar]^2\left|n,l,m_l\right&gt; \\
\hat{L}_z\left|n,l,m_l\right&gt; &amp;= L_z\left|n,l,m_l\right&gt; = m_l\hbar\left|n,l,m_l\right&gt;
\end{align}$$

Since no two orbitals have the same three quantum numbers, no two orbitals can have the same eigenvalues for all three operators. Consequently, there must be at least one operator which can be used to show that those two orbitals are orthogonal.