I'm given $\mathbb{R}^{n^2} = \mathbb{R}^{n \times n}$ with the usual metric $\langle A, B \rangle = \text{Trace} (AB^T) = \text{Trace} (A^T B)$.
Then $SO(n)$ is a submanifold of $\mathbb{R}^{n^2}$ of dimension $\frac{1}{2} n (n-1)$.
I have shown that the tangent space $T_{I} SO(n)$ consists of all skew-symmetric matrices with zero trace.
I wish to show that the normal space $N_I SO(n)$ consists of symmetric matrices. Furthermore, I need to calculate the second fundamental form tensor for $SO(n)$.
I'm given $X$ and $Y$ left-invariant vector fields determined by skew-symmetric matrices at the identity. If $O$ is an arbitrary orthogonal matrix, then $X(O) = OA $ and $Y(O) = OB$.
I need to prove that $$D_Y X(O) = \frac{1}{2} O [B, A] + \frac{1}{2} O (BA + AB) $$ and to find expressions for $\nabla_Y X$, $h(X,Y)$ and $[X, Y]$.
Here $\nabla_X Y$ denotes the covariant derivative on $SO(n)$, while $D_X Y$ denotes the covariant derivative on $\mathbb{R}^{n^2}$. Also, $h(X,Y)$ denotes the second fundamental form tensor.
Attempt:
I know that $D_X Y = \nabla_X Y + h(X,Y)$ is the Gauss formula.
On $\mathbb{R}^{n^2}$, I have $$ D_Y X (O) = D_{Y(O)} X = Y(O) X. $$ To calculate the latter, we take a curve through $O$ with tangent vector $Y(O) = OB$. Such a curve is $O \exp(Bt)$. Then $$ Y(O) X = \frac{d}{dt}_{t = 0} X (O \exp (Bt)) = OBA. $$
How can I prove that $$ 2 h(X,Y) = D_X Y + D_Y X $$ ? This would give me that $$h(X,Y) (O) = \frac{1}{2} O (AB + BA). $$
Also, I'm not sure how to calculate $\nabla_Y X$.
It seems to me that you've done roughly half of the work, and the part you're missing is seeing why $\nabla_XY|_O = \frac{1}{2}O[X,Y]$. This actually is a general result about Lie groups equipped with a bi-invariant metric (a metric which is both left- and right-invariant). A nice result that's often cited is that compact Lie groups are guaranteed a bi-invariant metric.
Proof: Letting $X,Y,Z \in \mathfrak{g}$ be left-invariant vector fields, we can compute $\left .\frac{d}{dt} \langle Ad_{\gamma(t)}Y, Ad_{\gamma(t)}Z\rangle \right |_{t=0}$, where $\gamma(t) = \exp(tX)$ and $Ad_{\gamma(t)} = \left (L_{\gamma(t)}\right )_*\circ \left (R_{\gamma(-t)}\right)_*$. Notice here that the bi-invariance of the metric yields
$$ \left . \frac{d}{dt} \langle Ad_{\gamma(t)}Y, Ad_{\gamma(t)}Z \rangle \right |_{t=0} \;\; =\;\; \left . \frac{d}{dt} \langle Y,Z\rangle \right |_{t=0} \;\; = \;\; 0 $$
where $\langle Y,Z\rangle$ is constant due to these being left-invariant vector fields. Observe on the other hand that \begin{eqnarray*} \left . \frac{d}{dt} \langle Ad_{\gamma(t)}Y, Ad_{\gamma(t)}Z\rangle \right |_{t=0} & = & \langle XY-YX, Z\rangle + \langle Y, XZ - ZX \rangle \\ & = & \langle [X,Y], Z\rangle + \langle Y, [X,Z]\rangle \;\; = \;\; 0. \end{eqnarray*} $$\tag*{$\Box$}$$
Proof: A general result from Riemannian geometry is that a Riemannian metric $\langle \cdot, \cdot\rangle$ satisfies the equation $$ \langle\nabla_XY,Z\rangle \;\; =\;\; \frac{1}{2} \left (X\langle Y,Z\rangle + Y\langle Z,X\rangle - Z\langle X,Y\rangle - \langle Y, [X,Z]\rangle - \langle Z,[Y,X]\rangle + \langle X, [Z,Y]\rangle \right ). $$
For each term of the form $X\langle Y,Z\rangle$, the action of $X$ is that of taking a derivative of the inner product $\langle Y,Z\rangle$ at each point. In the case of a bi-invariant Lie group with left-invariant vector fields $X,Y,Z \in \mathfrak{g}$, these terms all vanish, hence our equation in our special case reduces to
$$ \langle \nabla_XY, Z\rangle \;\; =\;\; \frac{1}{2} \left (-\langle Y, [X,Z]\rangle - \langle Z, [Y,X]\rangle - \langle X, [Z,Y]\rangle \right ). $$
By the first proposition we have the first and third terms both vanish, and we're left with
$$ \langle \nabla_XY, Z\rangle \;\; =\;\; -\frac{1}{2}\langle Z, [Y,X]\rangle \;\; =\;\; \left \langle Z, \frac{1}{2}[X,Y]\right \rangle \;\; =\;\; \left \langle \frac{1}{2}[X,Y], Z\right \rangle $$
where we used the anti-symmetry of the bracket and the symmetry of the metric. Since this equation is true for all left-invariant vector fields we conclude that $\nabla_XY = \frac{1}{2}[X,Y]$ on the left-invariant vector fields of $G$.
$$\tag*{$\Box$}$$
If we combine these two results with your computation that $D_XY(O) = OAB$, then we take the second-fundamental form $D_XY = \nabla_XY + h(X,Y)$. Noting the skew-symmetry of the bracket we have that \begin{eqnarray*} D_XY + D_YX & = & \nabla_XY + \nabla_YX + 2h(X,Y) \\ & = & \frac{1}{2}O[X,Y] + \frac{1}{2}O[Y,X] + 2h(X,Y) \\ & = & \frac{1}{2}O[X,Y] - \frac{1}{2}O[X,Y] + 2h(X,Y) \\ & = & 2h(X,Y). \end{eqnarray*}
From here you can conclude that $$ h(X,Y) \;\; =\;\; \frac{1}{2}O(XY + YX). $$
Note: A subtle point here is the fact that we assume $h(X,Y) = h(Y,X)$ even though in general the second fundamental form isn't necessarily a symmetric tensor. The reason why this works here is precisely because $N_ISO(n)$ consists of symmetric matrices, and by extension we have that $N_QSO(n) = \left \{QS \; | \; S = S^T\right \}$. This fact can actually be proven just with results about matrices and the trace. Since you've successfully proven that $T_ISO(n)$ is the space of skew-symmetric matrices, try showing that if $A \in T_ISO(n)$ then $Tr\left (A^TB\right ) = 0$ if and only if $B = B^T$.