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.
Best Answer
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.