Suppose I have a (simple, compact) Lie group $G$ with generators $T_i$. I would like to Taylor expand a (smooth) function $f:U\to\mathbb{R}$. I can always write a normal expansion like \begin{align} x &= e^{T_i \epsilon_i} x_0 \\ f(x) &= f(x_0) + \epsilon_i \left[\frac{\partial}{\partial \epsilon_i} f(e^{T_k\epsilon_k} x_0)\right]_{\epsilon=0} + \frac{1}{2}\epsilon_i\epsilon_j \left[\frac{\partial^2}{\partial \epsilon_i\epsilon_j}f(e^{T_k\epsilon_k}x_0)\right]_{\epsilon=0} + O(\epsilon^3) \end{align}
But now I would like to write this expansion in terms of the Lie derivative \begin{align} \nabla_i f(x_0) = \lim_{h \to 0} \frac{1}{h}\left(f(e^{hT_i}x_0)-f(x_0)\right) \end{align}
To first order, its simple: \begin{align} f(x) &= f(x_0)+e_i\nabla_i f(x_0) + O(\epsilon^2) \end{align}
But already the second order is not trivial, because \begin{align} \left[\frac{\partial^2}{\partial \epsilon_i\epsilon_j}f(e^{T_k\epsilon_k}x_0)\right]_{\epsilon=0} \neq \nabla_i\nabla_j f(x_0) \end{align} (The left expression is symmetric in $i,j$, the right is not, because Lie-derivatives do not commute, but instead have a relation like $\nabla_i\nabla_j-\nabla_j\nabla_i=-c_{ijk}\nabla_k$ with some structure constants $c$).
So what is the correct expression for higher terms of the taylor expansions using higher Lie-derivatives?