I struggle to reproduce a calculation from the Appendix of the paper "Anharmonic Oscillator: A Study of Perturbation Theory in Large Order", Physical Review D, 7 (6) 1973, [link to abstract][1]. 
I reproduce below my attempt. While main concern is the final result, I highlight two other doubts I have on their calculation in bold below. 

I need to compute the (first order approximation of the) probability current
$$J(x) = \frac{1}{2i}(-\Phi^{*}(x)\frac{d}{dx}\Phi(x)+\Phi(x)\frac{d}{dx}\Phi^{*}(x)$$
where $\bullet^{*}$ stands for the complex conjugate, for the WKB wave function
$$\Phi(x) \sim C_1 (x^2 - \epsilon x^4 - 4K - 2)^{-1/4} \times \exp \Big( -\frac{1}{2} \int _{x_0} ^{x} (t^2 - \epsilon t^4 -4K -2)^{1/2} dt \Big) $$
where $\epsilon$ is a (small) constant, $K$ the quantum number, $x_0 = (4K +2)^{1/2}$.
Let us define also
$x_1 \sim \frac{1}{\sqrt{\epsilon}}$ the distant turning point of the double well potential $U(x) = x^2 - \epsilon x^4$

Following the paper, they proceed as follows, 

> "We would now like to allow $x$, the upper endpoint of the WKB integral, to be larger than $x_1$. But of course, this expression for the WKB wave function is no longer valid as we pass the turning point at $x_1$. To avoid this difficulty we approach the point $x_1$ along a path which goes around $x_1$. <...> Fortunately, the $\textit{integral}$ in the WKB wave function depends only on the endpoint. Thus, for simplicity it may be taken entirely along the real axis. We then break the integral into two parts: (a) the portion below $x_1$ which is real <...> (b) the portion above $x_1$ which is imaginary"

Incidentally, I am still to understand their statement: **I do not get why the approach they advocate is legitimate, as it seems to me to defy the need for the connection formulas**.
However, this aspect put aside on the basis their argument is correct, I still cannot reproduce their calculation for $J$. 
This is my attempt.

Breaking the WKB integral as suggested
$$\Phi(x) \sim C_1 (x^2 - \epsilon x^4 - 4K - 2)^{-1/4} \times \exp \Big( -\frac{1}{2} \int _{x_0} ^{x_1} (t^2 - \epsilon t^4 -4K -2)^{1/2} dt \Big) \exp \Big( -\frac{1}{2} \int _{x_1} ^{x} (t^2 - \epsilon t^4 -4K -2)^{1/2} dt \Big) $$
For convenience, define

$$ E_0 ^1 = \exp \Big( -\frac{1}{2} \int _{x_0} ^{x_1} (t^2 - \epsilon t^4 -4K -2)^{1/2} dt \Big) $$ which is a real constant, and
$$ E(x) = \exp \Big( -\frac{1}{2} \int _{x_1} ^{x} (t^2 - \epsilon t^4 -4K -2)^{1/2} dt \Big)$$
$$ V(x) = U(x) -4K -2$$

allowing to re-write the WKB function as
$$\Phi(x) \sim C_1 V(x)^{-1/4}  E_0 ^1 E(x) $$
its derivative as
$$ \frac{d}{dx}\Phi(x) \sim  C_1 E_0 ^1 \Big( \frac{d}{dx} (V(x)^{-1/4}) E(x) + V(x)^{-1/4} \frac{d}{dx}E(x) \Big)$$
Now, 
$$ \frac{d}{dx} V(x)^{-1/4} = (-1/4) V(x)^{-5/4}\frac{d}{dx} V(x) $$
$$ \frac{d}{dx}E(x) = -\frac{1}{2} (x^2 - \epsilon x^4 -4K -2)^{1/2} $$
As $V(x)$ is negative for $ x > x_1$, it seems to make sense to neglect the second term in the bracket of the expression for $\frac{d}{dx}\Phi(x)$ getting to
$$ \frac{d}{dx}\Phi(x) \sim  C_1 E_0 ^1 \Big( \frac{d}{dx} \Big(V(x)^{-1/4}\Big) E(x)  \Big)$$
In the paper, **they claim they look for a leading order approximation of $J$, neglecting terms $ O(x^{-2}) < O (\epsilon)$ compared to $O(1)$, which is another statement I do not understand**: $\epsilon$ is a constant, what order should have?.
Anyhow, the complex conjugates are
$$\Phi^{*}(x) \sim C_1 E_0 ^1 \Big( V(x)^{-1/4}\Big)^{*}   E(x)^{*} $$
$$ \frac{d}{dx}\Phi(x)^{*} \sim  C_1 E_0 ^1 \Big( \frac{d}{dx} \Big(V(x)^{-1/4}\Big)^{*} E(x)^{*}  \Big)$$

Putting all in the expression for $J$, I get

$$ J(x) = \frac{1}{2i} \Big[-C_1^2 (E_0 ^1)^2  \Big( V(x)^{-1/4}\Big)^{*}  E(x)^{*}    \Big( \frac{d}{dx} \Big(V(x)^{-1/4}\Big) E(x) + C_1^2 (E_0 ^1)^2 V(x)^{-1/4}  E(x) \Big( \frac{d}{dx} \Big(V(x)^{-1/4}\Big)^{*} E(x)^{*}  \Big) \Big]$$

Now, if I understand correctly $E(x)$ is imaginary, so $ E(x)\Big(E(x) \Big)^{*} = 1$, and I get

$$ J(x) = \frac{1}{2i} \Big[-C_1^2 (E_0 ^1)^2  \Big( V(x)^{-1/4}\Big)^{*}      \Big( \frac{d}{dx} \Big(V(x)^{-1/4}\Big) + C_1^2 (E_0 ^1)^2 V(x)^{-1/4}   \Big( \frac{d}{dx} \Big(V(x)^{-1/4}\Big)^{*}   \Big) \Big]$$ or better rre-grouped

$$ J(x) = \frac{1}{2i} \Big\{-C_1^2 (E_0 ^1)^2 \Big [ V(x)^{-1/4}    \frac{d}{dx} \Big(V(x)^{-1/4}\Big)^{*} -   \Big( V(x)^{-1/4}\Big)^{*}     \frac{d}{dx} \Big(V(x)^{-1/4}\Big)  \Big]   \Big\}$$

The fact is, the result they report is simply

$$ J(x) = \frac{1}{2} C_1^2 (E_0 ^1)^2 = \frac{1}{2} C_1^2 \exp \Big( - \int_{x_0}^{x_1}(t^2 - \epsilon t^4 -4K -2)^{1/2} dt\Big) $$

I cannot figure it out, where is that I go wrong, any help would be greatly appreciated.


  [1]: https://doi.org/10.1103/PhysRevD.7.1620