Note that this is more general than the usual orthogonal matrix has orthogonal eigenvectors or symmetric matrix has orthogonal eigenvectors in that $A$ need not be orthogonal or symmetric, just square. Also this proof requires both directions ($\Leftrightarrow$).
I tried using the fact that a symmetric matrix has orthogonal eigenvectors by showing that ${A^T}A$ is a symmetric matrix and then showing that the eigenvectors of the original matrix A ($v_1$ and $v_2$) or some multiplication thereof with a matrix are eigenvectors of ${A^T}A$. To show the reverse direction is then easy.
We know that an eigenvalue of a square matrix is also an eigenvalue of its transpose and that an eigenvalue of a product of square matrices (A.B) is also an eigenvalue of the reverse product (B.A). The eigenvalues of real symmetric matrices are all real. Maybe this is not the right path, feel free to suggest another method.
This statement is mentioned page 37 of Gilbert Strang's book, Linear Algebra and Learning from Data.
After setting up the complete answer: it is also different from this in that the proof requires both directions (and here we only consider the real matrix $A$).
Decided to write up a detailed answer, that does not rely on any special theorems, and only relies on basic linear algebra properties, and some simple properties of the characteristic polynomial. The goal is to prove that $A$ is unitarily diagonalizale iff $A^*A = AA^*$. I've not seen many proofs like this, so here it is.
Let $A^*A = AA^*$ for a complex matrix $A\in M_n$. Such a matrix is called normal, and the properties discussed here generalize naturally to real normal matrices, in which case we replace the conjugate transpose $^*$ by $^T$. Then:
Proof: $\|Ax\|^2 = x^*A^*Ax = x^*AA^*x = \|A^*x\|^2$.
Proof: If $x\in N(A)$, then $0 = \|Ax\| = \|A^*x\|\implies x\in N(A^*)$.
Proof: If $y \in R(A), x\in N(A)$, then $x^*y = x^*A\xi = 0$, since $x\in N(A^*)$. $R(A)=R(A^*)$ is the corollary.
Proof: $(A + \lambda I)^*(A + \lambda I) = A^*A + \overline{\lambda} A+\lambda A^* + |\lambda|^2 I = (A + \lambda I)(A + \lambda I)^*$
Proof: $N(A - \lambda I) = N(A^* - \overline{\lambda} I)$ from 2. and 4.
Proof: Let $U = \begin{pmatrix}U_1 & U_2\end{pmatrix}$ be unitary, such that the columns of $U_1$ form a basis for $R(A)$ and the columns of $U_2$ form a basis for $N(A)$. Then using 2. and 3. and a calculation it is easy to verify that $U^*AU = B\oplus 0_{n-r}$. $B\in M_r$ implies that it is nonsingular and $$r = \rank(B) = \rank(B \oplus 0_{n-r}) = \rank(U^*AU) = \rank(A).$$
Proof: From 6. we know that $U^*(A - \lambda I)U = B\oplus 0_{n-r}$, thus $U^*AU = (B + \lambda I_r)\oplus \lambda I_{n-r} = C$, where $r = \rank(A - \lambda I)$. As such, the characteristic polynomial of $A$, $p_A(t) = p_C(t)$, the char. polynomial of $C$. But $p_C(t) = (\lambda - t)^{n-r}p_B(t - \lambda)$, and since $B$ is nonsingular, it follows that $\nu = n-r=\dim N(A - \lambda I)$.
Proof: Let $x,y$ be eigenvectors of $A$, corresponding to $\lambda,\mu$, respectively. Then using 5. $\lambda x^*y = x^*Ay = \mu x^*y \implies (\lambda - \mu)x^*y = 0\implies x^*y = 0$.
Proof: If $\nu_i$ is the algebraic multiplicity of $\lambda_i$, then we know that $n = \nu_1+\cdots+\nu_k$. This, combined with 7. and 8. gives the final result.
Now for the $\implies$ direction, if $A$ has $n$ orthonormal eigenvectors, then $A^*A = AA^*$.
Proof: Let $U\in M_n$ have the eigenvectors of $A$ as its columns, then it is unitary, $A = U\operatorname{diag}(\lambda_1,\ldots,\lambda_n)U^* = U\Lambda U^*$, and thus $$A^*A = U\overline{\Lambda} U^*U\Lambda U^* = U\Lambda \overline{\Lambda} U^* = AA^*.$$