How to understand the Second covariant derivative?

I am reading P. Li's lectures on Geometric analysis. On page 14, the author defines the second covariant derivative as follows:

Let $f$ be a smooth function on $M$. $\omega_1, \cdots, \omega_n$ be a local orthonormal basis of $T^*M$ around a fixed point $p$. And $d\omega_i=\sum_j\omega_{ij}\wedge\omega_j$. Then $$df=\sum_i f_i\omega_i$$ THen the author give the definition: $$f_{ij}\omega_j=df_i+f_j\omega_{ji}.$$

My question is how to understand this definition compare with the usual one, i.e. $$\nabla^2_{X,Y}f=XY(f)-\nabla_XY (f)$$

Also the third order covaiant derivative defines in the similar manner: $$f_{ijk}\omega_k=df_{ij}+f_{kj}\omega_{ki}+f_{ik}\omega_{kj}$$ The similar question for this expression, how to understand these indices?

Thank you for any detailed answer, I am kind of afraid of this local calculations compare with the ususal global definition. However, in most papers, the local calculations are more common. So I want to figure it out clearly.

They are actually the same. I am used to a slightly different notations and I will stick to mine. Let $\{\omega^i\}_{i=1}^n$ be a local orthonormal coframe dual to the orthonormal frame $\lbrace e_i\rbrace_{i=1}^n$ and we define the connection $1$-form $\omega_i^j$ by $\nabla _X e_i = \omega_i^j(X)e_j$ (Einstein notation). Then by differentiating $\delta_i^j=\omega^j(e_i)$, we have $\nabla _X \omega^j =-\omega_i^j (X)\omega^i$.
The equation $d\omega^i=\omega^j\wedge \omega_j^i$ (note that $\omega_i^j$ is $\omega_{ij}$ in Peter Li's note, and also $\omega_i^j=-\omega_j^i$ due to orthogonality), known as the (first) Cartan's structure equation, can be obtained by $$d\omega^i(X, Y)= (\nabla _X \omega^i)(Y)- (\nabla _Y \omega^i)(X)=-\omega^i_j(X)\omega^j(Y)+\omega^i_j(Y)\omega^j(X)=(\omega^j\wedge \omega ^i_j)(X, Y).$$ Alternatively, from the above it's easy to see that if $\omega^i_j$ satisfies the Cartan's equation, then $\nabla _X \omega^j =-\omega_i^j (X)\omega^i$. The connection 1-forms are exactly the analogue of the Christoffel symbols in a coordinates system.
Recall that for a function $f$, the covariant derivative $\nabla f$ coincides with $df$, i.e. $\nabla f(X)= df(X)(=X(f))$. The second covariant derivative (i.e. Hessian) is given by $\nabla ^2 f :=\nabla (\nabla f)= \nabla (df)$. In an orthonormal frame, we have \begin{equation*} \begin{split} f_{ij}:=(\nabla ^2 f)(e_i , e_j)= (\nabla (\nabla f))(e_i, e_j):= (\nabla _{e_i}\nabla f)(e_j)= &\nabla _{e_i}(\nabla _{e_j}f)- (\nabla f)(\nabla _{e_i}e_j)\\ =& e_i(e_j(f))-(\nabla _{e_i}e_j)f\\ =& e_i(f_j)-(\nabla _{e_i}e_j)f. \end{split} \end{equation*} where $f_j:= e_j(f)$. Using the language of moving frames, this computation becomes \begin{equation*} \begin{split} f_{ij}= e_i(f_j)-f_k \omega_j^k(e_i). \end{split} \end{equation*} As $dh=e_i (h) \omega^i$ for a function $h$, this means that $$f_{ij}\omega^i= df_j- f_k \omega^k_j= df_ j + f_k \omega_k^j.$$
It turns out that it doesn't matter if we define $f_{ij}$ as $\nabla ^2 f(e_i, e_j)$ or $\nabla ^2 f (e_j, e_i)$. However, for the third derivative, the order of $i, j, k$ does matter. In Peter Li's note, $f_{ijk}$ is defined as $(\nabla _k (\nabla ^2 f))(e_j, e_i)$, which in my convention is $(\nabla (\nabla ^2 f))(e_k, e_j, e_i)= :\nabla ^3 f(e_k , e_j, e_i)$. With this in mind, we now compute $f_{ijk}$.
As $\nabla ^2 f= f_{lm}\omega^l \omega^m$, we have \begin{equation*} \begin{split} f_{ijk}= (\nabla _k (\nabla ^2f))(e_j, e_i) =& \nabla _k (\nabla ^2 f (e_j , e_i))- \nabla ^2 f (\nabla _k e_j , e_i)-\nabla ^2 f ( e_j, \nabla _k e_i)\\ =&e_k(f_{ij})- \omega_j^l(e_k)f_{li}-\omega_i^l (e_k)f_{jl}\\ =& (df_{ij}- f_{li}\omega^l_k- f_{jl}\omega^l_i)(e_k) \end{split} \end{equation*} Therefore $$f_{ijk}\omega^k=df_{ij}- f_{li}\omega^l_k- f_{jl}\omega^l_i=df_{ij}+f_{li}\omega^k_l+f_{jl}\omega^i_l.$$