Wick rotation with magnetic fields How does Wick rotation work with magnetic fields?
Let us take up single-particle $d$+1 QM. Then the Euclidean time path integral is given (in $\natural$ units) by $$ \langle x| \exp(-t H) y\rangle = \mathcal{N} \int_{\gamma:[0,t]\to\mathbb{R}^d:\gamma(0)=x\land\gamma(t)=y}\exp(-\int_0^t\|\dot{\gamma}\|^2+V\circ \gamma)\mathrm{d}\gamma$$ where the Hamiltonian is given by $ H = -\Delta + V$ with $V:\mathbb{R}^d\to\mathbb{R}$ some potential. $\mathcal{N}$ is a normalization constant.
The decisive point here  was that we changed the Lagrangian from $$ L = \|\dot{\gamma}\|^2 - V\circ \gamma $$ to $$ L_E = \|\dot{\gamma}\|^2 + V\circ \gamma $$ in order to work in Euclidean time. What is nice about this integral is that the main contribution to it comes from the classical solution to the equations of motion (with inverted potential, however). In a sense then, its main use is to apply a Laplace approximation on it rather than actually calculate it. I.e., $$ \langle x| \exp(-t H) y\rangle \approx \exp(-\int_0^t\|\dot{\gamma_c}\|^2+V\circ \gamma_c)(1+\dots) $$ where $\gamma_c$ is the solution to the equation $$ \ddot{\gamma} = +(\nabla V)\circ \gamma \qquad(\gamma:[0,t]\to\mathbb{R}^d : \gamma(0)=x\land\gamma(t)=y)\,.$$

My question is, what is one supposed to do in order to follow this procedure in the presence of magnetic fields?

Allow me to elaborate to make the discussion precise. We take a magnetic field strength $B>0$ and pick the vector potential $$A(x) = \frac{1}{2} B e_3 \wedge x\qquad(x\in\mathbb{R}^d)$$ so that now our Hamiltonian is $$ H = (P-A)^2 + V \,.$$ It is known that in the presence of magnetic field, the Lagrangian should be now $$ L = \|\dot{\gamma}\|^2 - V\circ \gamma + \dot{\gamma}\cdot\left(A\circ\gamma\right) $$ and furthermore the Feynman-Kac-Ito formula does apply also in the presence of magnetic fields to yield $$ \langle x| \exp(-t H) y\rangle = \mathcal{N} \int_{\gamma:[0,t]\to\mathbb{R}^d:\gamma(0)=x\land\gamma(t)=y}\exp(-\int_0^t\|\dot{\gamma}\|^2+V\circ \gamma + \mathrm{i}  \dot{\gamma}\cdot (A\circ \gamma))\mathrm{d}\gamma \,. $$ However, now it is not clear at all what the classical solution even means. The equation of motion in Euclidean time $t\mapsto-\mathrm{i} t$ seems to yield a complex-valued path $\gamma$ (which is not what is being integrated on, so does one need to deform the contour integration in path space to reach this critical point?) which decays or explodes with time (in absence of $V$) unlike the usual cyclotron circular motion solutions.
So, in this case, what is the appropriate semiclassical approximation for $$ \langle x| \exp(-t H) y\rangle  \approx \quad ?$$
 A: Before I get to the core of the question, here are some peripheral comments about the path integral's definition:

*

*If the euclidean path integral converges when $A=0$, then it still converges when $A\neq 0$. The reason is simple: When $A=0$, the integrand $\exp\left(-\int dt\ L\right)$ is strictly positive. The $A$-term does not change the magnitude of $\exp\left(-\int dt\ L\right)$, so the $A$-term does not jeopardize the path integral's convergence. (The question already implicitly acknowledges this.)


*Even though the path integral converges, the $A$-term still needs to be handled with care. If we define the path integral by discretizing time, the $A$ factor must averaged over the interval that defines the discretized time-derivative. Otherwise, we get an incorrect factor of $2$ in the final result (reference 1, chapters 4 and 5).
Now I'll get to the core of the question:

... what is the appropriate semiclassical approximation...?

In this context, the name semiclassical approximation refers to an approximation in which the path integral is dominated by certain "paths." We don't need the "paths" to make any sense physically, because we're only using them as a tool to evaluate the path integral. This is already clear in the $A=0$ case, where the potential is inverted in the euclidean lagrangian. In the $A\neq 0$ case, the "classical solutions" are complex-valued, and that's okay. What matters physically is the path integral as a whole, Wick-rotated back to Lorentzian signature.
When the lagrangian is complex-valued, the approximation is usually called the saddle-point approximation instead of the semiclassical approximation. As far as I know, calculations have not been published for the specific case described in the question (a nonrelativistic particle in a background magnetic field), but the saddle-point approximation with complex-valued lagrangians has been the subject of some recent research. I'll review the approach that is being explored in the literature.
Using notation that is relatively standard in physics, the path integral can be written
$$
\newcommand{\bfx}{\mathbf{x}}
\newcommand{\bfA}{\mathbf{A}}
\newcommand{\pl}{\partial}
\int [d\bfx]\ \exp\big(-S[\bfx]\big)
\tag{1}
$$
with the euclidean action
$$
S[\bfx]\equiv
\int dt\ 
\left(\frac{\dot\bfx^2(t)}{2}+V\big(\bfx(t)\big)
 +i\dot\bfx(t)\cdot \bfA\big(\bfx(t)\big)
\right).
\tag{2}
$$
An overhead dot denotes a derivative with respect to $t$. In classical physics, the particle's coordinates $\bfx(t)$ would be regarded as smooth functions of time $t$, but in the context of the path integral, we can think of $t$ as an index instead. The path integral (1) involves an enormous (formally infinite) number of integration variables $\bfx(t)$ indexed by $t$. The action (2) is a function of these variables.
As a warm-up, consider a simpler integral with only one integration variable $x$:
$$
 \int dx\ \exp\big(-S(x)\big),
\tag{3}
$$
where the "euclidean action" $S(x)$ is a complex-valued function of the single real variable $x$. If $S(x)$ were real-valued and positive, then we could approximate the integral by expanding $S(x)$ about the point $x_0$ that minimizes $S(x)$. When $S(x)$ is complex, we can do this instead (section II in reference 3):

*

*Analytically continue $S(x)$ to be a holomorphic function of a complex variable $x$.


*Find the points $x$ where $S(x)$ satisfies the Cauchy-Riemann equations
$$
 \pl_x S(x)=0
\tag{4}
$$
with
$$
 \pl_x S\equiv 
 \frac{1}{2}
 \left(
 \frac{\pl S}{\pl\text{Re$(x)$}}
 -i\frac{\pl S}{\pl\text{Im$(x)$}}
 \right).
\tag{5}
$$
Points where equation (4) holds are called saddle points or just saddles, because the real part of $S(x)$ has a saddle-shape in the neighborhood of such a point. This is nicely illustrated in figure 1 of reference 3.


*Deform the contour of integration away from the real axis to a path that crosses through the saddle points in a special direction, namely the direction along which the real part of $S(x)$ increases most rapidly away from the saddle point. This is called a Lefschetz thimble. (Like much of the literature, reference 3 writes the integrand as $e^{S}$ instead of as $e^{-S}$, so the Lefschetz thimble is depicted as curving downward instead of upward.)
The result is a convergent approximation to the original integral (3). The direction along which the integration contour passes through the saddle points is important, because this optimizes the integral's convergence. Section 3 in reference 2 uses a simple example to illustrate the use of different contours.
This technique has been around for a long time, but apparently it has only recently been applied to path integrals, where the number of integration variables is enormous. Examples include references 4,5,6,7,8,9,10. I've only barely begun learning about the subject, but here are some general comments about how to implement it in the special case (1)-(2):

*

*The path integral (1) can be defined by discretizing time, deferring the continuous-time limit until the end of the calculation. After $t$ is discretized, it is manifestly only an index, and the number of integration variables $\bfx(t)$ is manifestly finite (but enormous). The saddle-point condition is analogous to (4), but now we require that condition with respect to all of the variables $\bfx(t)$ instead of just for one variable $x$. The integration "contour" is a submanifold of a many-dimensional complex space, chosen to pass through the saddle points in the direction that is most favorable for convergence.


*This requires analytically continuing the euclidean action $S[\bfx]$ to be a function of complex variables instead of real variables. If that seems too abstract, consider the case where the scalar potential $V(\bfx)$ and the vector potential $\bfA(\bfx)$ are both polynomials in the components of $\bfx$. Then the analytic continuation is relatively easy: simply replace each of the real variables by a complex variable.


*Even if we take $V(\bfx)$ and $\bfA(\bfx)$ to be polynomials, this still includes some interesting cases, such as the quartic potential $V(\bfx)\propto (\bfx\cdot\bfx)^2$ and the constant magnetic field $\bfA\propto (-x_2,\, x_1,\,0)$. Taking them to be polynomials is not required, but it makes the analytic continuation step easy, so I'd recommend trying it with polynomials first.
I haven't done the calculation myself, so I can't comment on the result. As a substitute, here is one interesting lesson from an explicit calculation in reference 11. The imaginary part of the euclidean action doesn't affect the magnitude of the integrand in (1), so we might be tempted to approximate the path integral by using the minima of the real part of the euclidean action. After all, the imaginary part can only improve the path integral's convergence. However, reference 11 uses an exactly-solvable example to demonstrate that this naive approach can give wrong answers. Deforming the integration "contour" (submanifold) to pass through the complex saddle points along Lefschetz thimbles is important.
By the way, terms that retain a factor of $i$ in the euclidean action are common in quantum field theory. Examples include:

*

*The time-derivative term in path integral for Dirac fermions. When relating the path integral formulation to the canonical formulation, this term can/should be regarded as part of the definition of the inner product, so we should be relieved that its coefficient is not affected by Wick rotation.


*The so-called theta term $\int F\wedge F$ in QCD in four-dimensional spacetime. Section 6.1 in reference 12 mentions an interesting connection (an analogy, at least) between theta terms and complex saddles.


*The Chern-Simons term in odd-dimensional spacetime. The fact that its coefficient is not affected by Wick rotation is important in the study of chiral anomalies, because the anomaly (which is matched by the Chern-Simons term through the anomaly-inflow mechanism) can only occur in the phase of the partition function, even when the signature is euclidean.
Applications of the complex-saddle-point method transcend these situations, though: the idea of using approximations based on complex saddles can be useful even when the euclidean action is not already complex.


*

*Schulman (1981), Techniques and Applications of Path Integration (Wiley-Interscience)


*Witten, Analytic Continuation Of Chern-Simons Theory (https://arxiv.org/abs/1001.2933)


*Feldbrugge, Lehners,and Turok, Lorentzian Quantum Cosmology (https://arxiv.org/abs/1703.02076)


*Witten, A New Look At The Path Integral Of Quantum Mechanics (https://arxiv.org/abs/1009.6032)


*Harlow, Maltz, and Witten, Analytic Continuation of Liouville Theory (https://arxiv.org/abs/1108.4417)


*Tanizaki, Nishimura, and Kashiwa, Evading the sign problem in the mean-field approximation through Lefschetz-thimble path integral (https://arxiv.org/abs/1504.02979)


*Behtash, Dunne, Schäfer1, Sulejmanpasic, and Ünsal, Complexified path integrals, exact saddles and supersymmetry (https://arxiv.org/abs/1510.00978)


*Cristoforetti, Di Renzo, Mukherjee, and Scorzato, Quantum field theories on the Lefschetz thimble (https://arxiv.org/abs/1312.1052)


*Cherman, Dorigoni, and Ünsal, Decoding perturbation theory using resurgence: Stokes phenomena, new saddle points and Lefschetz thimbles (https://arxiv.org/abs/1403.1277)


*Dunne and Ünsal, What is QFT? Resurgent trans-series, Lefschetz thimbles, and new exact saddles (https://arxiv.org/abs/1511.05977)


*Alexanian, MacKenzie, Paranjape, and Ruel, Path integration and perturbation theory with complex Euclidean actions (https://arxiv.org/abs/0802.0354), and also Problems With Complex Actions (https://arxiv.org/abs/hep-th/0609146)


*Dunne and Ünsal, New Methods in QFT and QCD (https://arxiv.org/abs/1601.03414)
A: *

*In this answer, we first of all want to point out that under a Wick rotation $t_E=it_M$ to Euclidean spacetime it is not just the magnetic part of the Lorentz force that turns imaginary. The electric part does as well. This is partly due to that the electric potential $A_0^M=iA_0^E$ is also Wick-rotated, cf. e.g. this Phys.SE post.


*In detail, the Lagrangian in Minkowski ($M$) spacetime is
$$\begin{align} L_M~=~&  \frac{m}{2}\left(\frac{d\vec{x}}{dt_M}\right)^2-U, \cr U ~=~&V -q \left(A^M_0 + \vec{A}\cdot \frac{d\vec{x}}{dt_M}\right),\end{align}\tag{A} $$
while the Lagrangian in Euclidean ($E$) spacetime reads
$$ \begin{align} L_E~=~&  \frac{m}{2}\left(\frac{d\vec{x}}{dt_E}\right)^2 +U, \cr U ~=~& V- iq \left(A^E_0 + \vec{A}\cdot \frac{d\vec{x}}{dt_E}\right).\end{align} \tag{B} $$
In contrast the scalar potential $V$ does not Wick rotate.


*Classical solutions $\vec{x}_0$ to Euler-Lagrange equations/Newton's 2nd law, i.e. stationary paths $\vec{x}_0$ for the action, which are real $\vec{x}_0\in \mathbb{R}^3$ in Minkowski signature, become complex $\vec{x}_0\in \mathbb{C}^3$ in Euclidean signature.


*We may then use complex function theory, Picard-Lefschetz theory & Lefschetz thimbles to deform the integration contour in the Euclidean path integral in order to only integrate over quantum fluctuations $\vec{\eta}$ along the direction of steepest descent
$$ \begin{align} \langle \vec{x}_f, t^E_f | \vec{x}_i, t^E_i \rangle 
~=~& \int_{\vec{x}(t^E_i)=\vec{x}_i}^{\vec{x}(t^E_f)=\vec{x}_f} \!  {\cal D}\frac{\vec{x}}{\sqrt{\hbar}} ~e^{-\frac{1}{\hbar}S_E[\vec{x}]}\cr
~=~& \int_{\vec{\eta}(t^E_i)=0}^{\vec{\eta}(t^E_f)=0} \! {\cal D}\vec{\eta} ~e^{-\frac{1}{\hbar}S_E[\vec{x}_0+\sqrt{\hbar}\vec{\eta}]},  \end{align}\tag{C} $$
and in this way achieve a WKB/semiclassical approximation around a complex classical path $\vec{x}_0$.
