Compute probability current from WKB approximation

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.
I reproduce below my attempt. While my 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}Big[-Phi^{*}(x)frac{d}{dx}Phi(x)+Phi(x)frac{d}{dx}Phi^{*}(x)Big]$$
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 and of unitary length, 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 re-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} dtBig) $$

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