Consider a (real) smooth manifold $M$, a (real, finite rank) vector bundle $E$ on $M$ and a connection $\nabla:\Gamma(TM)\times\Gamma(E)\to\Gamma(E)$. Let $X$ be a vector field on $M$. For simplicity I will assume that for some $t_0>0$ we have a flow $\theta:(-t_0,t_0)\times M\to M$ with infinitesimal generator $X$ although I believe the below can also be worded in terms of local flows. For a point $p\in M$ we have an integral curve of $\gamma_p:(-t_0,t_0)\to M$ of $X$ given by $\gamma_p(t)=\theta_t(p)$. We may consider the parallel transport map along this curve $P_t:E_p\to E_{\theta_t(p)}$.
Now, my first point is that it is easy to turn $P_t$ into a linear operator $\mathcal P_t$ on all of $\Gamma(E)$. Indeed, consider $\sigma\in\Gamma(E)$ and choose $\mathcal P_t(\sigma)$ so that for any $p\in M$ in the above notations we have $\mathcal P_t(\sigma)_{\theta_t(p)}=P_t(\sigma_p)$. Evidently, $\mathcal P_{-t}=\mathcal P_t^{-1}$. Furthermore, the identity $$\left.\frac d{dt}P_t^{-1}(\sigma_{\theta_t(p)})\right|_{t=0}=\nabla_X(\sigma)_p$$ can be now rewritten as $\frac d{dt}\mathcal P_{-t}|_{t=0}=\nabla_X$. The latter is an identity between linear operators on $\Gamma(E)$ and to compute the derivative one considers the pointwise (i.e. section-wise) convergence topology on $\mathrm{End}(\Gamma(E))$ induced by the pointwise convergence topology on $\Gamma(E)$.
Question. Is there anything inherently wrong with the above? If not, where can one find this approach of viewing parallel transport as an operator on the space of sections?
My second point is that in these terms one may obtain the formula for the Riemannian curvature operator. To me the below computation is a more convincing (or at least more precise) version of the usual hand-waving argument concerning the holonomy around an "infinitesimal parallelogram".
Choose vector fields $X$ and $Y$ both satisfying the above assumption and let $\mathcal P_t$ and $\mathcal Q_t$ be the corresponding operators on $\Gamma(E)$. Then, if we denote $\mathcal P_t=1+A(t)$ and $\mathcal Q_t=1+B(t)$, it is not hard to check that $$\left.\frac{\partial^2}{\partial t\partial s}\mathcal P_t^{-1}\mathcal Q_s^{-1}\mathcal P_t\mathcal Q_s\right|_{s=t=0}=\dot A(0)\dot B(0)-\dot B(0)\dot A(0)=\nabla_X\nabla_Y-\nabla_Y\nabla_X.$$ Note that this holds regardless of whether $X$ and $Y$ commute, i.e. whether the integral curves form a closed loop!
Finally, if we assume that $[X,Y]$ also satisfies the assumption and denote the corresponding operators by $\mathcal R_t$, we get $$\left.\frac{\partial^2}{\partial t\partial s}\mathcal R_{ts}\mathcal P_t^{-1}\mathcal Q_s^{-1}\mathcal P_t\mathcal Q_s\right|_{s=t=0}=\nabla_X\nabla_Y-\nabla_Y\nabla_X-\nabla_{[X,Y]}$$ so precisely the Riemannian curvature operator. The additional leftmost factor corresponds to the intuition of needing to correct by moving along $[X,Y]$ for time $ts$ to try to close the loop.
Question. Is there a mistake in this? If not, are such computations written out somewhere?
Update. Since this question has seen no activity, I reasked it on MathOverflow.