  In this paper, we present theoretical results on the statistical properties
of stationary, homogeneous and isotropic turbulence in incompressible flows in
three dimensions. Within the framework of the Non-Perturbative Renormalization
Group, we derive a closed renormalization flow equation for a generic $n$-point
correlation (and response) function for large wave-numbers with respect to the
inverse integral scale. The closure is obtained from a controlled expansion and
relies on extended symmetries of the Navier-Stokes field theory. It yields the
exact leading behavior of the flow equation at large wave-numbers $|\vec p_i|$,
and for arbitrary time differences $t_i$ in the stationary state. Furthermore,
we obtain the form of the general solution of the corresponding fixed point
equation, which yields the analytical form of the leading wave-number and time
dependence of $n$-point correlation functions, for large wave-numbers and both
for small $t_i$ and in the limit $t_i\to \infty$. At small $t_i$, the leading
contribution at large wave-number is logarithmically equivalent to $-\alpha
(\epsilon L)^{2/3}|\sum t_i \vec p_i|^2$, where $\alpha$ is a nonuniversal
constant, $L$ the integral scale and $\varepsilon$ the mean energy injection
rate. For the 2-point function, the $(t p)^2$ dependence is known to originate
from the sweeping effect. The derived formula embodies the generalization of
the effect of sweeping to $n-$point correlation functions. At large wave-number
and large $t_i$, we show that the $t_i^2$ dependence in the leading order
contribution crosses over to a $|t_i|$ dependence. The expression of the
correlation functions in this regime was not derived before, even for the
2-point function. Both predictions can be tested in direct numerical
simulations and in experiments.
