I'm following the proof in Goldstein's Classical Mechanics (ch 4-3) which instructs us to look at the double sum $a_{kl}a_{ki}a'_{ij}$, where we're using the Einstein notation for sums, and $a'$ are elements of the inverse matrix.
We then start from two extra self-evident equations $\delta_{li} a'_{ij}=a'_{lj}$ and similair $\delta_{kj} a_{kl}=a_{jl}$. We are supposed to arrive to $a_{jl}=a'_{lj}$ from these. However, the closest I can get to is $1/a_{jl}=a'_{lj}$ !
What I do is write down $a_{kl}\delta_{kj}=a_{jl}$, which is self evident, and substitute into it $a_{ki}a'_{ij}=\delta_{kj}$, which follows from the condition for the inverse.
Typically we define an orthogonal matrix $a_{ij}$ as one whose rows/columns form an orthonormal set of vectors under the dot product. You can write this condition as: $$ a_{ij}a_{ik} = \delta_{kj} $$ (with Einstein summation, this is the dot product of the $j^\text{th}$ and $k^\text{th}$ columns of the matrix). We can re-write the above as: $$ a_{ij}a^T_{ki} = [a^Ta]_{kj} = \delta_{kj} $$ Which (since matrix inverses are unique) identifies $a^T = a^{-1}$.