Time averages and differentiability Let $\varphi_t : M \rightarrow M$ be a smooth flow on a smooth manifold $M$. We may assume (although I'm not sure if this is important) that the flow preserves a smooth volume form on $M$. Given a continuous function $f : M \rightarrow \mathbb{R}$, is it true that if ALL the (partial) time averages
$$
A_T f(x) := {1 \over T} \int_0^T f(\varphi_t (x)) \, dt \hskip 1cm (T > 0)
$$
are smooth functions, then $f$ is also smooth? 
Remark. 
This is obviously true for the trivial action and it is non-trivially true for the case where $M = \mathbb{R}$ and the flow is the usual flow by translations. That's all I know for now. 
Motivation.
In the first volume of Dunford and Schwartz, where they motivate the Ergodic theorem (page 657), they state:
What is significant and measurable in the laboratory is not the quantity $f(\phi_t(x))$ but its average value
$$
 {1 \over T} \int_0^T f(\varphi_t (x)) \, dt
$$
computed over a certain time interval $0 \leq t \leq T$.
I'm asking whether knowing the regularity of all those average values says something about the regularity of the function, or "observable" $f$. It would be interesting to consider $f$ to be just locally integrable and then the question ties in a (very) little with Wiener's differentiation theorem. 
 A: $\textbf{Update with some corrections}$:
Unfortunately, my original answer had a serious gap. Sorry for the mix-up. The argument I gave only works when we have estimates on $A_T f(x)$ that are uniform in $T$. However, all the hard work was not useless as it made clear the properties that a counterexample would have. Without further ado, here is an example for which $A_T f(x)$ is smooth while $f(x)$ is not.
Let $M= \mathbb{R^2} =(t,x)$ with the flow $\phi_s (t, x) := (t+s, x)$.
Let $\Psi$ be the following function
$$\Psi(x) = \begin{cases} 
      0 & x\leq 0 \\
      e^{\frac{-1}{x^2}} &  x > 0 \\ 
\end{cases}$$
In particular, we want $\Psi$ to be a smooth function which vanishes to all orders at $0$. Let $f$ be the following function:
$$f(t,x)= \begin{cases}
 e^{\frac{-t^2}{(\Psi(x))^2}} & x >0 \\
 0 & x \leq 0
\end{cases} $$
The function $f$ is not smooth. In fact, it is not even continuous at $(0,0)$. To show that $A_T f(x)$ is smooth for any $T$, we must only show that it depends smoothly on $x$ near $(0,0)$. For fixed $T$>0, we have the following:
$$A_T(f(t_1,x))= \frac{1}{T} \int_{t_1}^{t_1+T} f(t,x) dt <  \frac{\sqrt{ \pi}}{T} \Psi(x) = o(x^k) ~ \forall k >0 $$
From this, it follows that $\frac{\partial^k A_T f}{\partial x^k}(t,0) \equiv 0$, so the derivative of $A_T f$ exists to all orders at $(0,0)$ for $T>0$. The smoothness of $A_T f(x)$ elsewhere is immediate because $f$ is smooth away from the origin.
For an example where $f(x)$ is continuous but not smooth, we can consider the function $\hat f(t,x) = x \cdot f(t,x)$. Here, the same analysis holds except that $f$ is continuous but not differentiable at $(0,0)$.
$\textbf{Original answer with the some corrections}$
This is not a full answer, but I believe that your colleague's result does all the heavy lifting away from the fixed points of the flow, at least once you have good estimates on $A_T f(x)$. The case near a fixed point of a flow seems more difficult. I'm not sure how that works exactly, but if one can understand it in the $1$-dimensional case (i.e. when the flow has a limit point), this might give some good insight. I broke the answer into several sections to help make it more readable.
$\textbf{Finding the right coordinates}$
In this case, the smooth flow locally foliates your manifold, so around any point $p$ one can use the inverse function theorem to choose local coordinates $\{ x_0, x_1, \ldots, x_{n-1} \}$ for which the flow acts as $\phi_t(x_0) = x_0+t$ and $\phi_t(x_i) = x_i$ otherwise. By Professor Burnol's argument (which I'm black-boxing for this argument) $f(  x_0, x_1, \ldots, x_{n-1})$ depends smoothly on $x_0$ when all the other coordinates are fixed.
What remains to show is that it depends smoothly on $x_1, \ldots, x_n$. We start with the case where $x_0$ is some fixed value, say 0. In this case, the $A_T f(x)$ is exactly the following:
$$ \frac{1}{T}\int_0^T f(t, x_1, \ldots, x_n) ~dt$$
For the rest of this, I will use the parameter $t$ for $x_0$ and refer to the rest of the coordinates collectively as $\mathbf{x}$.
We know that $A_T f(x)$ depends smoothly on $\mathbf{x}$ and further that $f$ depends smoothly on $t$. What remains to show is that $f$ depends smoothly on $\mathbf{x}$. We can't take derivatives with respect to $\mathbf{x}$, so we need to estimate some difference quotients. I will be using the Euclidean norm on the $\mathbf{x}$ coordinates for my estimates.
$\textbf{The uniform Lipschitz estimate}$
Take two points $\mathbf{x}^1$ and $\mathbf{x}^2$. Since $A_T f(x)$ is smooth, it is locally Lipschitz, and so we have the following estimate:
$$ \frac{1}{T \| \mathbf{x}^1- \mathbf{x}^2 \| } | \int_0^T f(t, \mathbf{x}^1) - f(t, \mathbf{x}^2)  ~dt | < \mathcal{C} $$
In order for the rest of the argument to work, we must assume that this estimate is uniform in the choice of $T$, $\mathbf{x}^1$ and $\mathbf{x}^2$. This allows us to pick $T>0$ small enough so that that for our fixed $\mathbf{x}^1$ and $\mathbf{x}^2$, we have the following:
$$ \| f(t, \mathbf{x}^i) -f(s, \mathbf{x}^i)  \|_{C^2} < \epsilon$$
for $0<t,s < T$.
The choice of $T$ may depend on $\mathbf{x}^1$ and $\mathbf{x}^2$. In particular, we don't expect to be able to do this uniformly in $\mathbf{x}$, but that's okay because we are only doing it for two points.
Then, we have the following estimate.
\begin{eqnarray}
| \frac{1}{T} \int_0^T \frac{ f(t, \mathbf{x}^1) - f(t, \mathbf{x}^2)}{\| \mathbf{x}^1- \mathbf{x}^2 \| }  ~dt | & \geq & | \frac{1}{T} 
\frac{ T(f(0,\mathbf{x}^1) -  f(0,\mathbf{x}^2))}{\| \mathbf{x}^1- \mathbf{x}^2 \|} | \\
& & - \frac{1}{T} \left( \frac{\partial f}{\partial t} (0, \mathbf{x}^1)+ \frac{\partial f}{\partial t} (0, \mathbf{x}^2) + 2\epsilon \right) \frac{T^2}{2} \\
& \geq & | \frac{f(0,\mathbf{x}^1) -  f(0,\mathbf{x}^2)}{\| \mathbf{x}^1- \mathbf{x}^2 \|}| - \left( \frac{\partial f}{\partial t} (0, \mathbf{x}^1)+ \frac{\partial f}{\partial t} (0, \mathbf{x}^2) + 2\epsilon \right) T \\
\end{eqnarray}
Letting $T$ be go to zero using the uniformity of the earlier estimate while still fixing $\mathbf{x}^1$ and $\mathbf{x}^2$ for now, our work shows the following.
$$ | \frac{f(0,\mathbf{x}^1) -  f(0,\mathbf{x}^2)}{\| \mathbf{x}^1- \mathbf{x}^2 \|}| \leq \mathcal{C} $$
Notice that this $\mathcal{C}$ is the same constant as before in the estimate of $A_T f(x)$. Repeating this for arbitrary pairs of $\mathbf{x}$, this gives a uniform Lipschitz estimate on $f$ (at $t=0)$ in terms of the $\mathbf{x}$ coordinates. 
$\textbf{A sketch of how to use induction for higher order regularity}$
In fact, we can repeat essentially the same argument but with three points $\mathbf{x}^1$, $\mathbf{x}^2$ and $ \lambda \mathbf{x}^1 + (1- \lambda) \mathbf{x}^2$ where instead of using the difference quotient, we use the second order difference quotient. This should yield a uniform $C^{1,1}$ estimate on $f$ in terms of the smoothness of $A_T f(x)$. With a uniform $C^{1,1}$ estimate, this shows that the function $f$ was actually differentiable in terms of $\mathbf{x}$ all along. You can induct on all orders to get smoothness of $f$ in $\mathbf{x}$. This shows that $f$ is smooth in both $t$ and $\mathbf{x}$.
