Let's try this. I didn't check all the estimates but the idea should be roughly as follows: since $M$ is compact, the flow $\Phi$ of $X$ is complete giving a one-parameter group action of $\mathbb{R}$ on $C^\infty(M)$. The latter is treated as Frechet space in the usual way. For a function $f \in C^\infty(M)$ define \begin{equation} f_s = \int_{\mathbb{R}} \Phi^*_\tau f e^{-1/\mu (\tau-s)^2} d\tau \end{equation} i.e. the convolution with the Gaussian. The integral yields again a smooth function since the derivatives (say covariant w.r.t. some auxilliary metric...) of $\Phi_t^*f$ are all bounded ($M$ is compact).

1. $f_0$ converges to $f$ for $\mu \to 0$ in the Frechet topology of $C^\infty(M)$. This is the standard kind of convolution argument. This shows that the span of the functions $f_s$ for varing $\mu$ and $s$ is dense.

2. $f_s$ has an entire extension in $s$. This should be OK from the explicit formula with the Gaussian. Replacing $s$ by $z \in \mathbb{C}$ gives a convergent integral since the factor $e^{-1/\mu t^2}$ dominates the other factors.

3. For $t \in \mathbb{R}$ one has $\Phi_t^*f_s = f_{t+s}$. This is just a change of variables in the integral.

4. This last point allows you to treat $z \mapsto f_z$ as an entire function with values in $C^\infty(M)$ extending the action of the flow on $f_0$.

Now, I guess this is as far as one can get. A few comments are perhaps necessary:

The ideas are sort of common with people doing Lie group representation theory. There, things are more complicated as there is typically nothing like a nice Gaussian available. This is very special to the abelian group $\mathbb{R}$.

The big question is of course what this is all good for? Well, I have no idea what your intention is to complexify the action. It is clear that there is no immediate geometric interpretation. On the other hand, the above construction uses only compactness but not the Riemannian metric, nice...

Let's try this. I didn't check all the estimates but the idea should be roughly as follows: since $M$ is compact, the flow $\Phi$ of $X$ is complete giving a one-parameter group action of $\mathbb{R}$ on $C^\infty(M)$. The latter is treated as Frechet space in the usual way. For a function $f \in C^\infty(M)$ define \begin{equation} f_s = \int_{\mathbb{R}} \Phi^*_\tau f e^{-1/\mu (\tau-s)^2} d\tau \end{equation} i.e. the convolution with the Gaussian.

Edit: Now the tricky point is whether the derivatives of $\Phi_t^*f$ stay bounded for all times $t$. Here the mere compactness will not help. The problem is that we will need something like that in order to differentiate into the integral to show that $f_s$ is still smooth.

Now the compactness of $M$ shows that the derivatives of the flow map on an interval $[0, 1]$ times $M$ we have bounded derivatives of $\Phi$. Using now the flow property $\Phi_n = \Phi_1 \circ \cdots \circ \Phi_1$ for all $n \in \mathbb{N}$ shows that the derivatives of $\Phi_t$ can grow at most exponentially in time, uniformly on $M$. Indeed, we can estimate the derivatives of $\Phi$ at $t$ by the derivatives of the $n$-fold composition of the map $\Phi_1$ and a sup over the derivatives of $\Phi$ on the above compact $[0, 1] \times M$ for some suitable $n$. This is a bit awkward to write down but should be OK (I hope). Now the integral together with the Gaussian can handle an exponential growth.

1. $f_0$ converges to $f$ for $\mu \to 0$ in the Frechet topology of $C^\infty(M)$. This is the standard kind of convolution argument. This shows that the span of the functions $f_s$ for varing $\mu$ and $s$ is dense.

2. $f_s$ has an entire extension in $s$. This should be OK from the explicit formula with the Gaussian. Replacing $s$ by $z \in \mathbb{C}$ gives a convergent integral since the factor $e^{-1/\mu t^2}$ dominates the other factors.

3. For $t \in \mathbb{R}$ one has $\Phi_t^*f_s = f_{t+s}$. This is just a change of variables in the integral.

4. This last point allows you to treat $z \mapsto f_z$ as an entire function with values in $C^\infty(M)$ extending the action of the flow on $f_0$.

Now, I guess this is as far as one can get. A few comments are perhaps necessary:

The ideas are sort of common with people doing Lie group representation theory. There, things are more complicated as there is typically nothing like a nice Gaussian available. This is very special to the abelian group $\mathbb{R}$.

The big question is of course what this is all good for? Well, I have no idea what your intention is to complexify the action. It is clear that there is no immediate geometric interpretation. On the other hand, the above construction uses only compactness but not the Riemannian metric, nice...