Exercise 95 (page 247) in Baez & Muniain's Gauge Fields, Knots and Gravity asks us to prove the well-known formula:
$$v - v' = \epsilon^2 F_{\mu \nu} v + o(\epsilon^2)$$
where $v \in E_p$ is a vector in the vector bundle $E \to M$ with a fixed connection, $v' \in E_p$ is the result of parallel transporting $v$ around a square of side $\epsilon$ in the $x^\mu$ - $x^\nu$ plane, and $F$ is the curvature tensor.
Specifically, the hint is to use the formula
$$ u(t) = \sum_{n=0}^{\infty} \left[ (-1)^n \int_{t \geq t_1 \geq \cdots \geq t_n \geq 0} A(\gamma'(t_1)) \cdots A(\gamma'(t_n)) dt_n \cdots dt_1 \right] u $$
which expresses the parallel transport of $u(0) = 0$ along the curve $\gamma(t)$ (it's possible to use the path-ordering operator to write this formula as an exponential, but I believe this more explicit version is the one that'll be useful for computation), keeping only terms up to order $\epsilon^2$.
Due to the linearity of holonomy, I believe it's sufficient to establish the desired formula for $v = e_i$ (with respect to some local trivialization that I'm supposing already fixed).
For convenience I recall
$$ F_{\mu \nu} e_i = \left( \partial_\mu A^j_{\nu i} - \partial_\nu A^j_{\mu i} + A^j_{\mu k} A^k_{\nu i} - A^j_{\nu k} A^k_{\mu i} \right) e_j $$
Let's WLOG suppose $p = 0$, and let $e^1_i$ be the result of parallel transporting $e_i$ from the origin to the first corner of the square (at $x^\mu = \epsilon, x^\nu = 0$). In this case, $\gamma' = \partial_\mu$. To first order in the paralell transport formula, we get
$$ e^1_i \approx e_i - \int_{0}^{\epsilon} A_\mu(0, \ldots, t_1, \ldots, 0) e_i dt_1$$
where $t_i$ belongs to the $\mu$ component. We could go to second order, but we'd get terms of the form $A_\mu A_\mu$ which seems a step in the wrong direction. To deal with this integral, I thought of approximating
$$A_\mu(0, \ldots, t_1, \ldots, 0) \approx A_\mu(0) + t_1 \partial_\mu A_\mu(0) $$
The first term integrates to give $-\epsilon A_\mu(0) e_i = -\epsilon A^j_{\mu i} e_j$, and the second gives $-(\epsilon^2/2) \partial_\mu A_\mu$ -- despite it being of second order, since it contains $\partial_\mu A_\mu$ it also doesn't look like it's going well. So if we unashamedly ignore that we'd get something like
$$e_i^1 \approx e_i - \epsilon A^j_{\mu i} e_j$$
The main problem is that by following this stratagem I don't see how we can possibly get the crossed derivative terms that need to show up in the final formula. For the 'horizontal' sides (along the $\mu$ direction) the integrals will always involve $A_\mu$ terms, and Taylor expanding those operators will give rise to $\mu$ derivatives of $A_\mu$ -- mutatis mutandis for $\nu$.
This is all to say I feel lost on how to go about this computation. Any help is greatly appreciated.