Reference for parallel transport around loop and its relation to curvature It is a well known fact that the geometric meaning of a linear connection's curvature can be realized as the measure of a change in a fiber element as it is parallel transported along a closed loop.
As far as I am aware, the formula goes as $$ F(X,Y)\sigma=-\lim_{t,s\rightarrow 0}\frac{1}{ts}\left(P^X_tP^Y_sP^X_{-t}P^Y_{-s}\sigma-\sigma\right), $$ where $P^X_t$ is parallel transport along $X$'s integral curves for time $t$, and I here assume that $[X,Y]=0$.
Unfortunately, none of my references contain a proof of this statement, at least not one that is useful to me now.
I can prove this statement using coordinate-based methods that often involve "expanding to first order" and other rather handwave-y methods (Weinberg, Wald etc. though Wald's procedure is actually fairly interesting and rigorous, it is not what I am looking for), but that's not what I am looking for.
I have also seen proofs of this statement (Lee's Manifolds and Differential Geometry) which involves starting with vectors $u$ and $v$ at points and extending them to vector fields in special ways, so that they are parallel along some curves, which simplifies things greatly. This is once again not what I am looking for, because it already presupposes that we know that $F$ is tensorial and that $F(X,Y)=\nabla_X\nabla_Y-\nabla_Y\nabla_X-\nabla_{[X,Y]}$.
In essence, what I am looking for is a derivation that $$ -\lim_{t,s\rightarrow 0}\frac{1}{ts}\left(P^X_tP^Y_sP^X_{-t}P^Y_{-s}\sigma-\sigma\right)=\nabla_X\nabla_Y\sigma-\nabla_Y\nabla_X\sigma $$ for commuting vector fields $X,Y$, without supposing we already know what curvature is (I need this for didactic reasons, for possible use in a general relativity lecture).
I have tried to do this on my own, but I am really terrible at formally manipulating these "flow" type of objects like $P^X_t$, and so far my attempts failed.
I would appreciate any reference, textbooks, lecture notes, articles, papers, which prove this formula (with the added caveat that I outlined in bold above). If it does so without assuming $X$ and $Y$ commute, it would especially be stellar, but I don't need that.
 A: A simple way is as follows.
Given a point $p$, the vector field
$(P_{-t}^X\sigma)(p)$ is obtained from $\sigma(\gamma(s))$ by backward parallel transport along the curve $\gamma(t):=\exp(tX)(p)$, up to time $0$. Denoting this backward parallel transport (from time $t$ to time $0$) by $Q_t^\gamma$, we have
$$ \frac{d}{dt}(P_{-t}^X\sigma)(p)=\frac{d}{dt}(Q_t^\gamma \sigma(t))=Q_t^\gamma(D_t\sigma(t))=(P_{-t}^X)(\nabla_X\sigma)(p), $$
where $\sigma(t):=\sigma(\gamma(t))$ is a vector field along $\gamma$ and $D_t$ is the covariant derivative along $\gamma$. So
$$P_{-t}^X\sigma=\sigma+\int_0^t P_{-t'}^X(\nabla_X\sigma)\,dt'=\sigma+t\nabla_X\sigma+O(t^2). $$
Hence,
$$\begin{align}&P_t^XP_s^YP_{-t}^XP_{-s}^Y\sigma=P_t^XP_s^YP_{-t}^X(\sigma+s\nabla_Y\sigma)+O(s^2)\\
&=P_t^XP_s^Y(\sigma+s\nabla_Y\sigma+t\nabla_X\sigma+ts\nabla_X\nabla_Y\sigma)+O(s^2)+O(t^2)\\
&=P_t^X(\sigma+t\nabla_X\sigma+ts\nabla_X\nabla_Y-ts\nabla_Y\nabla_X\sigma)+O(s^2)+O(t^2)\\
&=\sigma+ts\nabla_X\nabla_Y-ts\nabla_Y\nabla_X\sigma+O(s^2)+O(t^2)\end{align}$$
(of course the parallel transport applied to the error terms has the same order of magnitude, since it is obtained by solving a first-order ODE).
Now you just have to observe that the second order expansion of the smooth function $f(t,s):=(P_t^XP_s^YP_{-t}^XP_{-s}^Y\sigma)(p)$ has the form $\sigma(p)+tsV$ for some $V\in T_pM$: the coefficients of $s$, $t$, $s^2$ and $t^2$ all vanish since $f(t,0)\equiv\sigma(p)$ and $f(0,s)\equiv\sigma(p)$.
