Due to the relevance of the argument, after a nice discussion with a contribution of Carl Brannen, I decided to pursue this matter further. Indeed, the only way to have a covariant formulation of a gradient expansion is adding a time variable and taking the true time variable Wick rotated. In this way, for d=1+1 wave equation you will use d=2+1 wave equation and so on. In d=3+1 you will use d=4+1 wave equation. Let me explain with some equations what I mean. I consider again d=1+1 case as
but, instead to apply a gradient expansion to it, I apply this to the equation
being . As usual, I rescale time variable as and I take a solution series
Now I will get the set of equations
and so on. Let us note that, in this case, we can introduce two new spatial variables as and . These are conjugate variables as you know. So, already at the leading order I have solved my equation. Indeed, I note that
and so the Laplacian has the solution being f and g arbitrary functions. In this case the gradient expansion gives immediately the exact result making its application trivial as should be. Indeed, I take in the perturbation series and put and I get
that is the exact solution. Nice, it works! This means that a quantum field theory using gradient expansion exists and it is a strong coupling expansion. This result is surely less trivial than the one obtained above.