Solving Quantum Tunnelling Without Wick Rotation Edit
It seems that I haven't written my question clearly enough, so I will try to develop more using the example of quantum tunnelling. As a disclaimer, I want to state that my question is not about how to perform a Wick rotation in the path integral formulation!
Let's look a the probability of quantum tunnelling in the path integral formulation. The potential is given by $V[x(t)]=(x(t)^2-1)^2$, which has two minima at $x=\pm x_m=\pm1$. Given that the particle starts at $t=-\infty$ at $x=-x_m$, what is the probability that it is at $x=x_m$ at $t=\infty$. The  probability amplitude is given by 
$$K(x_m,-x_m,t)=\langle x_m|e^{-i \hat H t}|-x_m \rangle$$
The usual trick is to Wick rotate $t\to-i\tau$, compute everything in imaginary time using a saddle point approximation and at the end of the calculation rotate back to real time. I understand how it works. No problem with that.
What I want to understand is 


*

*how can I do the calculation without using the Wick rotation? 

*how does this solution connect to the Euclidean formulation?


In principle, we should be able to do the calculation with the path integral formulation in real time
$$\int Dx(t) e^{i S[x(t)]/\hbar}$$
In the stationary phase approximation we look for a complex path $x(t)$ which minimizes the action, and expand about this point.
Choose $m =1$ for simplicity. The equation of motion is
$$\ddot x-2 x+2x^3=0$$
which has no real solution, i.e. no Newtonian (classical) solution. But there is a complex function that solves it: $x_s(t)=i \,\tan(t)$. One problem is that it behaves pretty badly. If anyway I accept this a correct solution, I should be able to compute the gaussian fluctations, add up all the kinks/antikinks, etc. and recover the correct result (usually obtained with the euclidean action and  $\tau\to -it$). Am I right?
So my question is: is it possible to do the calculation that way, and if so, how is it related to the trick of going back and forth in imaginary time?
Original
I have a question on the mathematical meaning of the Wick rotation in path integrals, as it is use to compute, for instance, the probability of tunneling through a barrier (using instantons).
I am aware that when computing an ordinary integral using the Stationary Phase Approximation 
$\int dx e^{i S(x)/\hbar}$
with $x$ and $S$ real, one should look at the minimum of $S(z)$ in the whole complex plane, which can be for instance on the Imaginary axis.
In the case of a path integral, one wants to compute
$\int Dx(t) e^{i S[x(t)]/\hbar}$
and there is a priori no reason that the "classical path" from $x_a(t_a)$ to $x_b(t_b)$ (i.e. that minimizes $S[x(t)]$) should lie on the real axis. I have no problem with that. What I don't really get is the meaning of the Wick rotation $t\to -i\tau$ from a (layman) mathematical point of view, because it is not as if the function $x(t)$ is taken to be imaginary (say, $x(t)\to i x(t)$), but it is its  variable that we change !
In particular, if I discretize the path-integral (which is what one should do to make sense of it), I obtain
$\int \prod_n d x_n e^{i S(\{x_n\})/\hbar}$.
where $S(\{x_n\})=\Delta t\sum_n\Big\{ (\frac{x_{n+1}-x_n}{\Delta t})^2-V(x_n)\Big\}$
At this level, the Wick rotation applies on the time slice $\Delta t\to -i\Delta \tau$ and does not seem to be a meaningful change of variable in the integral
I understand that if I start with an evolution operator $e^{-\tau \hat H/\hbar}$ I will get the path integral after Wick rotation, but it seems to be a convoluted argument. 
The question is : Is it mathematically meaningful to do the Wick rotation directly at the level of the path-integral, and especially when it is discretized?
 A: I think the path-integral is a complete red herring here! I'll try to convince you that Wick rotation yields completely equivalent way of writing the Lagrangian in classical field theory.
Consider a classical action
$$S[x] = \int L[x(t)] dt$$
where $x:\mathbb{R} \to \mathcal{M}$ for some target manifold $\mathcal{M}$. The Lagrangian is schematically given by 
$$L[x(t)] = \left(\left.\frac{dx(s)}{ds}\right|_{s=t}\right)^2-V(x(t))$$
where $V(x(t))$ is some polynomial in $x(t)$ which (critically) involves no derivatives.
Now analytically continue $x$ to a function $\tilde{x}: \mathbb{C}\to\mathcal{M}$ by defining $\tilde{x}(c) = x(|c|)$ which is obviously analytic. Relabel $\tilde{x}$ as $x$ for simplicity. Define a new variable
$$\tau = it$$
inside the integral, and substitute. (Warning: there are mathematical subtleties about complex substitutions, which should be dealt with using Jordan's lemma). Ignoring the subtleties, the resulting integral is
$$S[x] = \int L[x(-i\tau)](-i)d\tau$$
Now let's examine what happens to the Lagrangian more carefully. Looking at the first term we have
$$\left(\left.\frac{dx(s)}{ds}\right|_{s=-i\tau}\right)^2$$
Change the differentiation variable to $u = is$ and this term becomes
$$\left(i\left.\frac{dx(u)}{du}\right|_{u=\tau}\right)^2$$
Relabelling $u\to s$ we see that the first term has the same form as originally, but with $t$ replaced by $\tau$ and an extra minus sign, viz.
$$-\left(\left.\frac{dx(s)}{ds}\right|_{s=\tau}\right)^2$$
Now on to the potential term. This is much simpler because $x(-i\tau)=x(|-i\tau|)=x(\tau)$ by definition so the potential term is just
$$V(x(\tau))$$
which is of exactly the same form as originally. Now we define a Euclidean Lagrangian
$$L_E(x(\tau)) = \left(\left.\frac{dx(s)}{ds}\right|_{s=\tau}\right)^2+V(x(\tau))$$
Putting it all together we find
$$S[x] = i\int L_E[x(\tau)]d\tau$$
Finally defining 
$$S_E[x] = \int L[x(t)]dt$$
we see that it's mathematically equivalent to calculate the path integral as
$$\int \mathcal{D}x\exp(iS[x]) \textrm{ or } \int \mathcal{D}x\exp(-S_E[x])$$
A: The usefulness of Wick rotation lies in the convergence properties of the path integral. If you look at the integrand of the path integral in Minkowski space, 
$$\int\mathcal{D}\phi\;e^{iS_M},$$
you can see that it is an oscillating function. The integral of an oscillating function can in general be considered problematic. Wick rotation, which is equivalent to a change from Minkowski to Euclidean space, changes this expression to
$$\int\mathcal{D}\phi\;e^{-S_{E}}.$$
The integrand is now an exponentially dampening function, which is well behaved compared to the original one. This enables one to do many calculations which would otherwise be hard to do. Furthermore, it is noteable that if one identifies this new "Euclidean time" with inverse temperature, the path integral corresponds to the partition function of statistical physics.
A: Others have already given correct answers, but perhaps I can help clear up your confusion by focusing on the simplest case, where your action is free and your discretization has just one point.
In this case, Wick rotation amounts to defining the integral $\int_{\mathbb{R}} e^{iax^2} dx$ to be the limit $\lim_{\alpha \to -ia} F(\alpha)$ where $F(\alpha) = \int_{\mathbb{R}} e^{-\alpha x^2} dx$.  
Note that there is no change of variable in the integral.  Instead you have a free parameter.  The integral is well-defined for real values of this parameter, and you define the integral at complex values via analytic continuation.
