Variation of a time-ordered exponential Consider the time-ordered exponential (Wilson line):
$$
U(t_{f},t_{i})
=
\mathcal{T}\text{exp}\left(-i\int_{t_{i}}^{t_{f}}\mathcal{A}(t)dt\right)\tag{1}
$$
Where $\mathcal{A}(t)$ is some matrix-valued function (gauge connection), and $\mathcal{T}$ denotes the time ordering.
I want to calculate the first order variation $\delta U(t_{f},t_{i})$, under a transformation:
$$
\mathcal{A}\mapsto\mathcal{A}+\delta\mathcal{A}\tag{2}
$$
There is probably a clever way to do this, but all I have been able to do so far is work with the definition. I have found the following.
$$\begin{align*}
U(t_{f},t_{i})\mapsto &
\sum_{n=0}^{\infty}
\frac{(-i)^{n}}{n!}\int_{t_{i}}^{t_{f}}\cdots \int_{t_{i}}^{t_{f}}\mathcal{T}\left((\mathcal{A}(t_{1}+\delta\mathcal{A}(t_{1})\cdots(\mathcal{A}(t_{n}+\delta\mathcal{A}(t_{n})\right)dt_{1}\cdots dt_{n}
\\
&=
\sum_{n=0}^{\infty}
\frac{(-i)^{n}}{n!}\int_{t_{i}}^{t_{f}}\cdots \int_{t_{i}}^{t_{f}}\mathcal{T}\left(\mathcal{A}(t_{1})\cdots\mathcal{A}(t_{n})+(\delta\mathcal{A})(t_{1})\mathcal{A}(t_{2})\cdots\mathcal{A}(t_{n})+\cdots\right)dt_{1}\cdots dt_{n}
\\
&=
U(t_{f},t_{i})
+
\sum_{n=0}^{\infty}
\frac{(-i)^{n}}{n!}\int_{t_{i}}^{t_{f}}\cdots \int_{t_{i}}^{t_{f}}\mathcal{T}\left((\delta\mathcal{A})(t_{1})\mathcal{A}(t_{2})\cdots\mathcal{A}(t_{n})+\cdots+\mathcal{A}(t_{1})\mathcal{A}(t_{2})\cdots\mathcal{A}(t_{n-1})(\delta\mathcal{A})(t_{n})\right)dt_{1}\cdots dt_{n}
\\ 
\end{align*}\tag{3}$$
Where we have dropped terms more than linear in the variations, and where each term has the form of a product of $\mathcal{A}(t_{i})$'s with a single term having been replaced by a variation $(\delta\mathcal{A})(t_{i})$.
The desired result is:
$$
\delta U(t_{f},t_{i})
=
-i\int_{t_{i}}^{t_{f}}U(t_{f},t)(\delta\mathcal{A})(t)U(t,t_{i})dt.\tag{4}
$$
I have been able to show that the two sides agree to second order, and I can kind of understand how moving the $(\delta\mathcal{A})(t_{j})$'s through the time ordered product will lead to $U(t_{f},t)$ appearing on the left and $U(t,t_{i})$ on the right, but after a few hours of messing around I am having trouble with the specifics. Any help would be much appreciated.
EDIT:
Using Qmechanic's answer I have been able to utilise the group property of the time-ordered exponentials and a discretisation of time to find:
$$
\delta U(t_{n},t_{1})=\sum_{i=1}^{n}U(t_n,t_{i+1})\delta U(t_{i+1},t_{i})U(t_{i},t_{1})
$$
Taking the interval $|t_{i+1}-t_{i}|<<1$ then lets us disregard time-ordering on this interval, so that:
$$
U(t_{i+1},t_{i})
\mapsto
\text{exp}\left(-i\int_{t_{i}}^{t_{i+1}}\mathcal{A}(t)+\delta\mathcal{A}(t)dt\right)
\approx
e^{(t_{i+1}-t_{i})(\mathcal{A}(t_{i})+\delta\mathcal{A}(t_{i}))}
$$
If we now assume that $\delta\mathcal{A}(t_{i})$ and $\mathcal{A}(t_{i})$ commute, it is easy to see that:
$$
\delta U(t_{i+1},t_{i})
=
-i(t_{i+1}-t_{i})U(t_{i+1},t_{i})\delta\mathcal{A}(t_{i})
$$
Which along with the group property, turns my sum into a Riemann sum and in the continuum limit gives the desired integral. However, I fail to see why this commutativity should hold. There is probably a pretty obvious general reason, but I can't seem to see it and would really appreciate some clarification.
 A: The time-ordered exponential
appears in the quantum mechanical time evolution operator
$$U(t_f,t_i)~=~T\exp\left[-\frac{i}{\hbar}\int_{t_f}^{t_i}\! dt~H(t)\right].\tag{A}$$
It satisfies a group property
$$U(t_3,t_1)~=~U(t_3,t_2)U(t_2,t_1).\tag{B}$$
Assume that we have a sufficiently fine discretization of time
$$t_n~>~t_{n-1}~>~\ldots~>~t_2~>~t_1,\tag{C}$$
so that
$$\forall t~\in~ [t_i,t_{i+1}]:~~ (t_{i+1}-t_i) || H(t) || ~\ll~ \hbar. \tag{D}$$
Next use the group property (B) and Leibniz rule to deduce that an infinitesimal variation is
$$\begin{align}\delta U(t_n,t_1)~=~&\sum_{i=1}^n U(t_n,t_{i+1})~\delta U(t_{i+1},t_i)~U(t_i,t_1)\cr
~\approx~&-\frac{i}{\hbar}\sum_{i=1}^n U(t_n,t_{i+1})~
(t_{i+1}-t_i) \delta H(t_i)~U(t_i,t_1)
\end{align} \tag{E}$$
In the continuum limit we heuristically get OP's desired formula
$$ \delta U(t_f,t_i)~=~-\frac{i}{\hbar}\int_{t_i}^{t_f}\! dt~U(t_f,t)~
 \delta H(t)~U(t,t_i)\tag{F}$$
For a similar formula without time-ordering, see this related Phys.SE post.
A: The time-ordered exponential can be written in the following form:
$$ U(t_f, t_i) = \lim_{N \to \infty} e^{-i A(t_1) \Delta t} \cdots e^{-i A(t_N) \Delta t} , $$
where $\Delta t = (t_f - t_i) / N$ and $t_k = t_i + k \Delta t$.
From there I find it quite intuitive that
$$ \begin{aligned} \delta U(t_f, t_i) &= \lim_{N \to \infty} \sum_{k=1}^N e^{-i A(t_1) \Delta t} \cdots e^{-i A(t_k) \Delta t} \bigl[-i\, \delta A(t_k)\, \Delta t\bigr]\, e^{-i A(t_{k+1}) \Delta t} \cdots e^{-i A(t_N) \Delta t} \\ &= \int_{t_i}^{t_f} d\tau\, U(\tau, t_i) \bigl[ -i\, \delta A(\tau) \bigr]\, U(t_f, \tau) \end{aligned} $$
by the definition of the Riemann integral.
