Theorem. Let $\textbf{P}$ be a permutation matrix corresponding to the permutation $\rho:\{1,2,\dots,n\}\to\{1,2,\dots,n\}$. Then $\textbf{P}^t=\textbf{P}^{-1}.$
Proof. First note the following identity for the kronecker delta:
$$\sum_k \delta_{ik}\delta_{kj} = \delta_{ij}$$
Now, by definition we have $\textbf{P}_{ij}=\delta_{i\rho(j)}$. Further, $\textbf{P}_{ij}^t = \delta_{\rho(i)j}$. But then:
$$(\textbf{PP}^t)_{ij}=\sum_k \textbf{P}_{ik}\textbf{P}^t_{kj} =\sum_k \delta_{i\rho(k)}\delta_{\rho(k)j} = \delta_{ij}$$
Hence $\textbf{PP}^t = I$ and we are done.$\,\square$
Is there a more elegant way to go about this proof (perhaps using $S_n$ somehow)? Is the $\delta$ definition of a permutation matrix for $\rho$ even the best definition? (I haven't in the past seen any others...)
Let $e_1,\dots,e_n$ be a fixed orthonormal basis of a real $n$ dimensional inner product space $V$.
Then $P$ can be regarded as the linear transformation $V\to V$ mapping $e_i\mapsto e_{\rho(i)}$.
Since $\rho$ is a permutation, we will have $$\langle Pe_i,Pe_j\rangle\ =\ \delta_{i,j}\ =\ \langle e_i,e_j\rangle$$ so, by linearity, $P$ will preserve the inner product of any two vectors. Finally, by the definition of adjoint (which, as matrix, corresponds to the transpose), the above gives $$\langle P^tPe_i,e_j\rangle=\langle e_i,e_j\rangle$$ for all $i,j$, and consequently, $P^tP$ must be the identity.