Calculating the Berry curvature in case of degenerate levels (Non abelian Berry curvature): issue The Berry phase accumulated on a path can be described by a matrix when we look at adiabatic time evolution with a Hamiltonian with degenerate energy levels.
The Berry phase matrix is given by 
$$
\gamma_{mn}= \int_\mathcal{C} \left\langle m(R) \right | \nabla_R \left| n(R) \right \rangle . d R.
$$
here $R$ parametrizes the said path and $ A_{mn}= \left\langle m(R) \right | \nabla_R \left| n(R) \right \rangle$. Now what I want to do is calculate the Berry Curvature, something that, if I assume my path above is closed and has three determining coordinates$R_1$, $R_2$ , $R_3$  is $\vec{F}$ such that 
$$
\mathbf{\gamma}=\int_\mathcal{S} \vec{F}.d\vec{s}
$$
note that the $\gamma$ and $F$ here are matrices and we're integrating over the surface $\mathcal{S}$ enclosed by curve $\mathcal{C}$ 
What is stopping me from applying the stokes theorem to $\gamma_{mn}$ and getting $\vec{F}_{mn}=\nabla_R\times A_{mn} $?
It is said that the answer contains a matrix commutator $[A_i,A_j]_{mn}$ c because this berry phase is non abelian. but I seem to be missing something fundamental.
Edit:
Note: this also corresponds to problem 2 Chapter 2 of Topological insulators and superconductors by Bernevig and Hughes
 A: The short answer is that only if the Berry curvature is defined by: (in matrix notation):
$$F_{\mu \nu} = \partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu} + [A_{\mu}, A_{\nu}]$$
it becomes gauge covariant, i.e., for a gauge transformation:
$$A_{\mu} \rightarrow g^{-1}A_{\mu} g+g^{-1}\partial_{\mu}g$$
$g \in U(N)$ ($N$ is the degeneracy of the level), the curvature transforms as:
$$F_{\mu \nu} \rightarrow g^{-1}F_{\mu \nu} g$$
The partial expression $F_{\mu \nu} = \partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ will be covariant (actually invariant) only with respect to Abelian gauge transformations and not the whole $U(N)$ group.
The reason that the Berry curvature cannot be obtained as in the Abelian case by a direct  application of the Stokes theorem is that the Stokes theorem does not exist(in the usual sense) in the non Abelian case becuse the holonomy needs to be path ordered:
$$\mathcal{Hol}(A) = \mathrm{P} e^{\int_C A_{\mu}dx^{\mu}}$$
However, it can be applied sequentially for an infinitesimal path. let us choose a square path in the $(xy)$ plane centered at $(x,y)$ with infinitesimal sides of length $\Delta x$ and $\Delta y$ as depicted in the picture.

Therefore
$$ \begin{matrix}
\mathcal{Hol}(A) \approx  e^{\int_a^b A_{\mu}dx^{\mu}} e^{\int_b^c A_{\mu}dx^{\mu}} e^{\int_c^d A_{\mu}dx^{\mu}} e^{\int_d^a A_{\mu}dx^{\mu}}
\\\approx (1+A_x\Delta x-\frac{1}{2}\partial_yA_x\Delta x \Delta y) (1+A_y\Delta y+\frac{1}{2}\partial_x A_y\Delta x \Delta y) (1-A_x\Delta x-\frac{1}{2}\partial_yA_x\Delta x \Delta y)(1-A_y\Delta y+\frac{1}{2}\partial_x A_y\Delta x \Delta y) 
\\ \approx 1+(\partial_x A_y - \partial_y A_x + [A_x, A_y])\Delta x \Delta y\approx (1+ F_{xy}\Delta x \Delta y)
\end{matrix} $$
Please notice that the commutator term is created due to the need to exchange the order of the second and the third expressions in the product in order to cancel $A_x \Delta_x$.
