How do we use contour integration to invert a Laplace transform?

Question

A handy way of inverting Laplace transforms is the use of the Bromwich integral combined with the use of the Cauchy residue theorem. Below I will demonstrate an example where I can invert the Laplace transform analytically obtaining a closed form expression. On the other hand by using the Cauchy theorem I will get a completely wrong result which confuses me a lot. I would like to understand where is the error in my reasoning.

---

Let us take $x > b >0 $, $\theta < 1$, $\mu_0 <1/2 $ and $\nu = (1-2 \mu_0)/(2(1-\theta)) $ and then let us consider a following function below. We have:

\begin{eqnarray} u_\alpha^{(b)}(x) &:=& E_x \left[ e^{-\alpha \tau_b} \right] \\ &=& \left( \frac{x}{b} \right)^{\frac{1}{2} - \mu_0} \cdot \frac{K_\nu\left[ \sqrt{2} \frac{x^{1-\theta}}{1-\theta} \sqrt{\alpha}\right]}{K_\nu\left[ \sqrt{2} \frac{b^{1-\theta}}{1-\theta} \sqrt{\alpha}\right]} \tag{1} \end{eqnarray}

Note that the quantity in $(1)$ is a valid Laplace transform. Indeed we know from here that for small values of $ z $ the modified Bessel function of the second kind behaves as $ K_\nu(z) \simeq z^{-\nu} $ for $ \nu > 0 $. As such it is clearly seen that $ \lim_{\alpha \rightarrow 0} u_\alpha^{(b)}(x) = 1 $ as it should be.

It turns out ---- Ito, K.; McKean, H. P. jun., Diffusion processes and their sample paths, Berlin-Heidelberg-New York: Springer-Verlag. XVII, 321 p. (1965). ZBL0127.09503. section 4.6 page 130. ---- that the quantity above is a Laplace transform of the distribution of the first hitting time $ \tau_b := inf(t>0 | X_t > b) $, given $X_0 = x$, in an one-dimensional diffusion process of the kind $d X_t = \mu_0 X_t^{2 \theta-1} dt + X_t^\theta d B_t $ where $\theta \in [1/2,1) $ with $B_t$ being the Brownian motion.

Let us now take $\nu = 1/2 $ , which implies $\mu_0 = \theta/2$, and let us compute the inverse Laplace transform in two different ways.

First way:

We define ${\mathfrak A}(x,b) := (x^{1-\theta} - b^{1-\theta})/(1-\theta) $ and we firstly compute the Bromwich integral analytically as below:

\begin{eqnarray} \rho_{\tau_b}(t) &=& \frac{1}{2\pi \imath} \int\limits_{\imath {\mathbb R}} e^{\alpha \cdot t} \cdot u_\alpha^{(b)}(x) d\alpha \\ &=& \frac{1}{2\pi} \left( \frac{b}{x} \right)^{\mu_0-\theta/2} \int\limits_{\mathbb R} e^{\imath \alpha t} e^{-\sqrt{2} \sqrt{\imath \alpha} {\mathfrak A}(x,b)} d \alpha \\ &=& \frac{1}{2\pi} \left( \frac{b}{x} \right)^{\mu_0-\theta/2} 2 Re\left[ \int\limits_0^\infty 2 \beta e^{\imath (\beta^2 t - \beta {\mathfrak A}(x,b))} e^{-\beta {\mathfrak A}(x,b)} d\beta \right] \\ &=& \left( \frac{b}{x} \right)^{\mu_0-\theta/2} \frac{1}{\sqrt{2 \pi}} t^{-3/2} {\mathfrak A}(x,b) e^{-\frac{{\mathfrak A}(x,b)^2}{2 t}} \end{eqnarray}

In[1952]:= {b, z} = RandomReal[{1/2, 3/2}, 2]; alpha =.;
x = RandomReal[{2 b, 10}];
th = RandomReal[{1/2, 1}];
nu = 1/2;
mu0 = 1/2 - nu (1 - th);
NIntegrate[
  Exp[I alpha z] (
   x^(1/2 - mu0)
     BesselK[nu, (Sqrt[2] Sqrt[(I alpha)] x^(1 - th))/(1 - th)])/(
   b^(1/2 - mu0)
     BesselK[nu, (Sqrt[2] Sqrt[(I alpha)] b^(1 - th))/(
     1 - th)]), {alpha, -Infinity, Infinity}]/(2 Pi)
(b/x)^(mu0 - th/2)
  NIntegrate[
   Exp[I alpha z] E^(((b x)^-th (Sqrt[2] Sqrt[
       I alpha] (b^th x - b x^th)))/(-1 + th)) , {alpha, -Infinity, 
    Infinity}]/(2 Pi)
(b/x)^(mu0 - th/
  2) 2 Re[NIntegrate[
    2 beta Exp[
      I beta^2 z] E^(-( ((1 + I) beta (x^(1 - th) -  b^(1 - th)))/(
      1 - th))) , {beta, 0, Infinity}]]/(2 Pi)
A =  ( (x^(1 - th) -  b^(1 - th)))/(1 - th);
(b/x)^(mu0 - th/
  2) 2 Re[NIntegrate[ 
    2 beta Exp[I (beta^2 z - beta A)] E^- ( beta A), {beta, 0, 
     Infinity}]]/(2 Pi)
(b/x)^(mu0 - th/2) (A E^(-(A^2/(2 z))) Sqrt[1/(2 \[Pi])])/z^(3/2)
During evaluation of In[1952]:= Infinity::indet: Indeterminate expression 0 [Infinity] encountered.
During evaluation of In[1952]:= Infinity::indet: Indeterminate expression 0 [Infinity] encountered.
Out[1957]= 0.00352682 + 2.03297*10^-13 I
Out[1958]= 0.00352682 + 2.03297*10^-13 I
Out[1959]= 0.00352682
Out[1961]= 0.00352682
Out[1962]= 0.00352682

In the first line from the top we just inserted the definition of the Laplace transform into the integral, in the second line we substituted $\beta = \sqrt{\alpha}$ and finally we did the integral analyticaly. It is easy to see that the result is a true probability distribution being supported on the positive half-axis. It integrates to one.

Second way:

We use the Cauchy residue theorem as applied to the contour below:

Bear in mind that we took $\nu = 1/2$. Then the modified Bessel function in question is given in closed form and we know that the Laplace transform $u_\alpha^{(b)}(x) $ does not have any poles in the negative complex half plane. Therefore we can just write:

\begin{eqnarray} \frac{1}{2\pi \imath} \int\limits_{\Gamma_1} e^{\alpha \cdot t} \cdot u_\alpha^{(b)}(x) d\alpha + \frac{1}{2\pi \imath} \int\limits_{\Gamma_2} e^{\alpha \cdot t} \cdot u_\alpha^{(b)}(x) d\alpha = 0 \\ \Leftrightarrow \rho_{\tau_b}(t) + \frac{1}{2\pi} \left( \frac{b}{x} \right)^{\mu_0-\theta/2} \int\limits_{\pi/2}^{3/2\pi} e^{R e^{\imath \phi} t} e^{-\sqrt{2} \sqrt{R e^{\imath \phi}} {\mathfrak A}(x,b)} \imath R e^{\imath \phi} d\phi = 0 \end{eqnarray} Now, clearly the second integral above disappears in the limit $R \rightarrow +\infty$ because $Re[ R e^{\imath \phi} ] <0 $ for $\phi \in [\pi/2,3/2\pi] $.

This would imply that our probability distribution is just identically zero which is obviously wrong.

Where is then the error in my reasoning?

Answer 1

Owing to to Steven Gubkin's comment above and also to a great blog post on those things we can give an answer straight away. Here we do this for $\nu= 1/2$ only but it might be possible to do it in general ; we simply need to identify all the branching points in the later case.

Let us assume that $\nu= 1/2$. We apply the Cauchy theorem to the contour below where the distance of the vertical line in the positive half-plane from the imaginary axis is $c$ and the radius of the circle is $R$. We have:

Again, by means of the Cauchy theorem the sum of the five integrals below is equal to zero. We have:

\begin{eqnarray} \int\limits_{c - \imath \sqrt{R^2-c^2}}^{c + \imath \sqrt{R^2-c^2}} d\alpha e^{-\sqrt{2} {\mathfrak A} \sqrt{\alpha}} \cdot e^{\alpha t} + \underline{\int\limits_{\frac{\pi}{2} - \arcsin(\frac{c}{R})}^\pi R\imath e^{\imath \phi} d\phi \cdot e^{-\sqrt{2} {\mathfrak A} \sqrt{R e^{\imath \phi}}} \cdot e^{R e^{\imath \phi} t}} + \int\limits_0^R e^{-\sqrt{2} {\mathfrak A} \imath \sqrt{\alpha}} \cdot e^{-\alpha t} d\alpha - \int\limits_0^R e^{+\sqrt{2} {\mathfrak A} \imath \sqrt{\alpha}} \cdot e^{-\alpha t} d\alpha + \underline{\int\limits_{\pi}^{\frac{3 \pi}{2} + \arcsin(\frac{c}{R})} R\imath e^{\imath \phi} d\phi \cdot e^{-\sqrt{2} {\mathfrak A} \sqrt{R e^{\imath \phi}}} \cdot e^{R e^{\imath \phi} t}} = 0 \end{eqnarray}

I think that everything above is pretty straightforward except for one thing, i.e. that as we go around the small circle next to the origin we gain a phase, i.e. $\sqrt{\alpha} \rightarrow - \sqrt{\alpha} $ and that is why the integrals along the horizontal intervals next to the negative real axis do not cancel each other.

Now let us analyze the second an the fifth terms. We have:

\begin{eqnarray} &&2nd + 5th= \\ &&\int\limits_{\frac{\pi}{2} - \arcsin(\frac{c}{R})}^{\frac{3\pi}{2} + \arcsin(\frac{c}{R})} R\imath e^{\imath \phi} d\phi \cdot e^{-\sqrt{2} {\mathfrak A} \sqrt{R e^{\imath \phi}}} \cdot e^{R e^{\imath \phi} t} \underbrace{=}_{\zeta = e^{\imath \pi}} \\ &&\int\limits_{\frac{c}{R} + \imath \sqrt{1- \frac{c^2}{R^2}}}^{\frac{c}{R} - \imath \sqrt{1- \frac{c^2}{R^2}}} d\zeta e^{-\sqrt{2}{\mathfrak A} \sqrt{R \zeta}} \cdot e^{R \zeta t} \underbrace{=}_{R\rightarrow \infty}\\ && \int\limits_{\frac{c}{R} + \imath}^{\frac{c}{R}-\imath} d \zeta e^{R \zeta t} = \\ && e^{c t} (-2 \imath) \frac{\sin(R t)}{R t} \underbrace{=}_{R \rightarrow \infty} \\ && (-2\imath \pi) e^{c t} \delta(t) \end{eqnarray}

We can see that except for $t=0$ the contribution from the arcs in the contour vanish. We are therefore left with three terms only. In the limit $R \rightarrow \infty$ we have:

\begin{eqnarray} &&\int\limits_{c - \imath \infty}^{c + \imath \infty} e^{-\sqrt{2} {\mathfrak A} \sqrt{\alpha}} \cdot e^{\alpha t} d\alpha = \int\limits_0^\infty \left[ e^{+\sqrt{2} {\mathfrak A} \imath \sqrt{\alpha}} - e^{-\sqrt{2} {\mathfrak A} \imath \sqrt{\alpha}}\right] e^{-\alpha t} d\alpha \\ % 1st &&= 2 \imath \cdot e^{- \frac{{\mathfrak A}^2}{2 t}} \cdot Im \left[ \int\limits_0^\infty e^{-t(\beta - \frac{{\mathfrak A}}{\sqrt{2} t} \imath)^2} 2 \beta d\beta\right] \\ % 2nd &&= 2 \imath \cdot e^{- \frac{{\mathfrak A}^2}{2 t}} \cdot Im \left[ \int\limits_{{\mathbb R}_+ } e^{-t \beta^2} 2\left[\beta + \frac{{\mathfrak A}}{\sqrt{2} t} \imath\right] d\beta + \int\limits_{\imath [-\frac{{\mathfrak A}}{\sqrt{2} t}, 0] } e^{-t \beta^2} 2\left[\beta + \frac{{\mathfrak A}}{\sqrt{2} t} \imath\right] d\beta \right] \\ %3 rd &&= 2 \imath \cdot e^{- \frac{{\mathfrak A}^2}{2 t}} \cdot Im\left[ \int\limits_{{\mathbb R}_+} e^{-t \beta^2} 2 \frac{{\mathfrak A}}{\sqrt{2} t} \imath d\beta \right] \\ &&= 2 \imath \cdot e^{- \frac{{\mathfrak A}^2}{2 t}} \cdot \frac{{\mathfrak A}}{\sqrt{2} t^{3/2}} \cdot 2 \underbrace{\int\limits_0^\infty e^{-\beta^2} d\beta}_{ \frac{\sqrt{\pi}}{2}} \\ && = 2 \imath \cdot e^{- \frac{{\mathfrak A}^2}{2 t}} \cdot \frac{{\mathfrak A}}{\sqrt{2} t^{3/2}} \cdot \sqrt{\pi} \\ && (2\pi \imath ) \cdot e^{- \frac{{\mathfrak A}^2}{2 t}} \cdot \frac{{\mathfrak A}}{t^{3/2}} \cdot \frac{1}{\sqrt{2\pi}} \end{eqnarray}

Note: The step from the second to the third line involves a substitution $\beta - {\mathfrak A} \imath/(\sqrt{2} t) \rightarrow \beta$ and then using the Cauchy theorem applied to an infinite strip in the forth quadrant.

Then we need to multiply the result by the prefactor $(b/x)^{\mu_0 - \theta/2} $ and we retrieve the pdf of the first hitting time from the body of the question. q.e.d.

Answer 2

Here we provide an answer for arbitrary $\nu > 0 $. We use the same contour as in my answer above. Now in order to proceed we need to identify the poles and the branching point of the Laplace transform $u_\alpha^{(b)}(x) $. Since the modified Bessel function of the second kind is multi-valued at the origin only there will be one branching point, namely at $\alpha= 0 $. In addition to this there might be poles, which are zeros of the denominator. From here we learn that in the right half-plane (positive real part) there are no zeros. All the zeros of the function in question are located in the interior of the second and the third quarters.This means that the potential poles of the Laplace transform have to satisfy:

\begin{equation} \sqrt{2} \frac{b^{1-\theta}}{1-\theta} \cdot \sqrt{\alpha} = z \tag{1} \end{equation} where $|Arg[z]| > \pi/2$ (with the convention that $Arg[z] \in [-\pi,\pi]$. But equation $(1)$ does not have any solutions in $\alpha$ (this is because the square root of a complex number has an argument that is always between $-\pi/2$ and $\pi/2$. Therefore there are no poles.

On the other hand the contribution from the branching point (at the origin) comes from two terms, the first one being an integral of the integrand over $(-R,0)$ and the second one being an integral of the integrand evaluated at $\sqrt{\alpha} \rightarrow - \sqrt{\alpha} $, an integral taken over $(0,-R)$.

This gives us the following formula below:

\begin{eqnarray} \rho_{\tau_b}(t) &=& &&\left( \frac{x}{b} \right)^{\nu(1-\theta)} \frac{1}{\pi} \int\limits_0^\infty Im\left[ \frac{K_\nu[\sqrt{2} \frac{x^{1-\theta}}{1-\theta} (-\imath) \sqrt{\alpha}] }{K_\nu[\sqrt{2} \frac{b^{1-\theta}}{1-\theta} (-\imath) \sqrt{\alpha}] } \right] \cdot e^{-\alpha t} d\alpha \tag{1} \end{eqnarray}

Below we fixed $x,b,\nu$(as in the caption of the plot) and we plotted the probability distribution of the first hitting time $\tau_b$ as a function of time for different values of $\theta=linspace(0.6,09,5)$ (from violet to red respectively). Here the numerical calculation (based on the Bromwich integral) and the analytical result (based on $(1)$) are on the left and on the right respectively. As we can see there is a perfect match .

(*Run (*Procedure for finding zeros *) before running the below.*)
(*For the Bromwich integral to exist we must have: x> b. This is \
because BesselK[z] ~ Exp[-z]/Sqrt[z] as |z| --> Infinity and |arg[z]| \
< 3/2 Pi.*)
(Now we need to satisfy the condition (6a), i.e. the Laplace 

transform below  must vanish at x-->-Infinity. This is equivalent to 

Cos[Pi (1-th)] >0  which is equivalent to Pi \in [1/2,3/2] - 2 Z.)
ths = Array[# &, 5, {6/10, 9/10}];
cols = Table[
   ColorData["Rainbow", i/(Length[ths] - 1)], {i, 0, Length[ths] - 1}];
(mu0 = 1/2-nu(1-th);)
Clear[LT]; 
LT[alpha_, b_, th_] := (
 x^(nu (1 - th))
   BesselK[nu, (Sqrt[2] Sqrt[(alpha)] x^(1 - th))/(1 - th)])/(
 b^(nu (1 - th))
   BesselK[nu, (Sqrt[2] Sqrt[(alpha)] b^(1 - th))/(1 - th)]);
zs = Array[# &, 50, {1/2, 2}];
(Compute numerically.)
vals = Table[
   NIntegrate[
       Exp[I alpha #] LT[I alpha, b, 
         ths[[ii]]], {alpha, -Infinity, +Infinity}]/(2 Pi) & /@ 
    zs, {ii, 1, Length[ths]}];
vals1 = Table[(x/b)^(nu (1 - ths[[ii]])) ( (1/Pi) NIntegrate[
         2 beta Im[ 
           BesselK[
            nu, (Sqrt[2]  x^(1 - ths[[ii]]))/(
             1 - ths[[ii]]) (-I) beta]/ 
           BesselK[
            nu, (Sqrt[2] b^(1 - ths[[ii]]))/(
             1 - ths[[ii]]) (-I) beta ]] Exp[-beta^2 #], {beta, 0, 
          Infinity}]) & /@ zs, {ii, 1, Length[ths]}];
str = ToString[{x, b, nu}];
b =.; pl1 = 
 ListPlot[Table[Transpose[{zs, vals[[ii]]}], {ii, 1, Length[ths]}], 
  PlotStyle :> cols, 
  AxesLabel :> {Subscript[[Tau], b], 
    Subscript[[Rho], Subscript[[Tau], b]]}, 
  PlotLabel :> "Bromwich integral", ImageSize :> 400];
pl2 = ListPlot[
   Table[Transpose[{zs, vals1[[ii]]}], {ii, 1, Length[ths]}], 
   PlotStyle :> cols, 
   AxesLabel :> {Subscript[[Tau], b], 
     Subscript[[Rho], Subscript[[Tau], b]]}, 
   PlotLabel :> "Cauchy theorem", ImageSize :> 400];
GraphicsGrid[{{pl1, pl2}}, PlotLabel -> "x,b,[Nu]=" <> str]

Update:

Formula $(1)$ isn't really a solution to the problem because it still involves an improper integral. However, we are going to argue that if $\nu = 1/2 + {\mathbb N} $ a closed form expression can be obtained. As an example let us take $\nu = 3/2$. Then the result reads as follows:

\begin{eqnarray} &&\rho_{\tau_b}(t) = \frac{{\mathfrak A}(x,b) e^{-\frac{{\mathfrak A}(x,b)^2}{2 t}} }{\sqrt{2 \pi } \sqrt{t}} \cdot \\ && \left( % \frac{(\frac{x}{b})^{1-\theta}}{t} + % \right. \\ && \left. (1-\theta)^2 b^{2 \theta-2} \left( -1 + % \sqrt{\frac{\pi }{2}} \sqrt{t} (1-\theta )^1 b^{\theta -1} e^{\frac{\left(t (\theta -1) b^{\theta }-b {\mathfrak A}(x,b)\right)^2}{2 b^2 t}} \text{erfc}\left(\sqrt{\frac{\left(b {\mathfrak A}(x,b)+t (1-\theta ) b^{\theta }\right)^2}{2 b^2 t}}\right) \right) % % \right) \tag{2} \end{eqnarray} where ${\mathfrak A}(x,b):= (x^{1-\theta} - b^{1-\theta})/(1-\theta)$.

(*nu=3/2*)
{b} = RandomReal[{1/2, 4}, 1];
x = RandomReal[{2 b, 10}];
th = RandomReal[{1/2, 1}];
AA[x_, b_] :=    (x^(1 - th) - b^(1 - th))/(1 - th);
ts = Array[# &, 100, {1/10, 10}];
nu = 3/2;
vals0 = (x/b)^(nu (1 - th)) ( (1/Pi) NIntegrate[
        2 beta Im[ 
          BesselK[nu, (Sqrt[2]  x^(1 - th))/(1 - th) (-I) beta]/ 
          BesselK[nu, (Sqrt[2] b^(1 - th))/(
            1 - th) (-I) beta ]] Exp[-beta^2 #], {beta, 0, 
         Infinity}]) & /@ ts;
vals0 = NIntegrate[(x/b)^((1 - th) nu) 1/
      Pi (2 beta Im[ 
         BesselK[nu, (Sqrt[2]  x^(1 - th))/(1 - th) (-I) beta]/ 
         BesselK[nu, (Sqrt[2] b^(1 - th))/(
           1 - th) (-I) beta ]] Exp[-beta^2 #]), {beta, 0, 
      Infinity}] & /@ ts;
vals1 =
  Exp[-(AA[x, b]^2/(2 #))] AA[x, b]/Sqrt[
     2 [Pi] #] (b^(-2 + th)   (-(1 - th)^2 b^th + (b x^(1 - th))/#) +
        b^(-3 + 3 th) Sqrt[#] (1 - th)^3 Sqrt[[Pi]/2]
         E^((b^th # (-1 + th) - b AA[x, b])^2/(2 b^2 #))
          Erfc[(b^th # (1 - th) + b AA[x, b])/(b Sqrt[2 #])]) & /@ ts;
Abs[vals0/vals1 - 1] // Mean
b =.; ListPlot[Transpose[{ts, #}] & /@ {vals0, vals1}, 
 PlotRange :> All, PlotStyle :> {{Purple, Diamond}, {Blue, Triangle}},
  AxesLabel :> {Subscript[[Tau], b], 
   Subscript[[Rho], Subscript[[Tau], b]]}]

Now, by using the asymptotic expansion for the complementary error function it is easily shown that :

\begin{eqnarray} \lim\limits_{t \rightarrow \infty}\rho_{\tau_b}(t) = \frac{{\mathfrak A}(x,b) e^{-\frac{{\mathfrak A}(x,b)^2}{2 t}} }{\sqrt{2 \pi } \sqrt{t}} \cdot \left( \frac{b^{-2 \theta } x^{-2 \theta } \left(b^{\theta +1} x^{\theta +1}+x^2 b^{2 \theta }+b^2 x^{2 \theta }\right)}{t^2 (\theta -1)^2} + O(\frac{1}{t^3}) \right) \end{eqnarray}

which implies that fractional moments exist if and only if the order of the moments is strictly smaller than three half, i.e. $E[ (\tau_b)^\zeta ] < \infty $ iff $ \zeta < 3/2 $.

Update 1:

Here we evaluate formula $(1)$ for $\nu = 1/2 + n$ where $n\in{\mathbb N}$.We define ${\tilde x}(\beta): = \sqrt{2} \frac{x^{1-\theta}}{1-\theta} \cdot \beta$ and ${\tilde b}(\beta): = \sqrt{2} \frac{b^{1-\theta}}{1-\theta} \cdot \beta$ and we have:

\begin{eqnarray} \rho_{\tau_b}(t) &=& \left( \frac{x}{b} \right)^{\nu(1-\theta)} \frac{1}{\pi} \int\limits_0^\infty Im\left[ \frac{K_\nu[ (-\imath){\tilde x}(\beta) ] }{K_\nu[ (-\imath) {\tilde b}(\beta)] } \right] \cdot 2 \beta e^{-\beta^2 t} d\beta \\ &=& \left( \frac{x}{b} \right)^{\nu(1-\theta)} \frac{1}{\pi} \int\limits_0^\infty \frac{J_\nu({\tilde b}(\beta)) \cdot Y_\nu({\tilde x}(\beta)) - J_\nu({\tilde x}(\beta)) \cdot Y_\nu({\tilde b}(\beta))}{J_\nu({\tilde b}(\beta))^2 + Y_\nu({\tilde b}(\beta))^2 } \cdot 2 \beta e^{-\beta^2 t} d\beta \\ &=& \left( \frac{x}{b} \right)^{\nu(1-\theta)} \frac{Re}{\pi} \int\limits_0^\infty \left( \frac{{\tilde b}(\beta)}{{\tilde x}(\beta)} \right)^{n+\frac{1}{2}} \cdot e^{\imath ({\tilde x}(\beta)-{\tilde b}(\beta))} \cdot \frac{ {\mathfrak F}_\nu({\tilde x}(\beta),{\tilde b}(\beta)) } { \sum\limits_{j=0}^n (n-j+1)^{(2 j)} \frac{(2j-1)!!}{2^j j!} \cdot {\tilde b}(\beta)^{2(n-j)} } \cdot 2 \beta e^{-\beta^2 t} d\beta \\ &=& \frac{Re}{\pi} \int\limits_0^\infty \frac{ {\mathfrak F}_\nu({\tilde x}(\beta),{\tilde b}(\beta)) } { \sum\limits_{j=0}^n (n-j+1)^{(2 j)} \frac{(2j-1)!!}{2^j j!} \cdot {\tilde b}(\beta)^{2(n-j)} } \cdot 2 \beta e^{-\beta^2 t + \imath \sqrt{2} {\mathfrak A}(x,b) \cdot \beta} d\beta \\ &=& e^{-\frac{{\mathfrak A}(x,b)^2}{2 t}} \frac{2}{\pi} Re \int\limits_0^\infty \frac{ {\mathfrak F}_\nu({\tilde x}(\beta),{\tilde b}(\beta))\cdot \beta } { \sum\limits_{j=0}^n (n-j+1)^{(2 j)} \frac{(2j-1)!!}{2^j j!} \cdot {\tilde b}(\beta)^{2(n-j)} } \cdot e^{-t \cdot\left( \beta - \frac{\imath {\mathfrak A}(x,b)}{\sqrt{2} t} \right)^2} d\beta \\ &=& e^{-\frac{{\mathfrak A}(x,b)^2}{2 t}} \frac{2}{\pi} Re \sum\limits_{j=0}^{n+1} \frac{d^j_\beta }{j!} \left.[{\mathfrak F}_\nu({\tilde x}(\beta),{\tilde b}(\beta))\cdot \beta]\right|_{\beta \rightarrow 0} \cdot \\ && \sum\limits_{\xi=1}^{2 n} \underbrace{ \frac{[\zeta_\xi^{(n)}]^j}{\prod\limits_{\eta=1,\xi \neq \eta }^{2 n} (\zeta_\xi^{(n)} - \zeta_\eta^{(n)})} }_{{\mathfrak C}^{(n,j)}_\xi} \cdot \int\limits_0^\infty \frac{e^{-t \cdot\left( \beta - \frac{\imath {\mathfrak A}(x,b)}{\sqrt{2} t} \right)^2}}{(\sqrt{2} \frac{b^{1-\theta}}{1-\theta} \cdot \beta - \zeta_\xi^{(n)})} d\beta \\ &=& e^{-\frac{{\mathfrak A}(x,b)^2}{2 t}} \frac{2}{\pi} Re \sum\limits_{j=0}^{n+1} \frac{d^j_\beta }{j!} \left.[{\mathfrak F}_\nu({\tilde x}(\beta),{\tilde b}(\beta))\cdot \beta]\right|_{\beta \rightarrow 0} \cdot \\ && \sum\limits_{\xi=1}^{n} {\mathfrak C}^{(n,j)}_\xi \cdot \underline{ \int\limits_0^\infty e^{-t \cdot\left( \beta - \frac{\imath {\mathfrak A}(x,b)}{\sqrt{2} t} \right)^2} \left[ \frac{1}{(\sqrt{2} \frac{b^{1-\theta}}{1-\theta} \cdot \beta - \zeta_\xi^{(n)})} + \frac{(-1)^{j+1}}{(\sqrt{2} \frac{b^{1-\theta}}{1-\theta} \cdot \beta - \bar{\zeta_\xi^{(n)}})} \right] d\beta } \\ \end{eqnarray} In the first line we substituted $\beta = \sqrt{\alpha}$ and in the second line we used the identity $K_\nu(-\imath x) = \pi/2 \imath^{1+\nu} (J_\nu(x) + \imath Y_\nu(x))$. In the third line we evaluated the fraction in the integrand in a following way. The denominator was being obtained by using a corresponding identity from here and the numerator was being obtained by using the respective spherical Bessel function identities here and here.

Now the ${\mathfrak F}_\nu(\cdot,\cdot)$ function is a certain multinomial in its arguments and are defined as follows:

\begin{eqnarray} {\mathfrak F}_\nu(x,b) := (-\imath) \sum\limits_{l_1=0}^n \sum\limits_{l_2=0}^n {\mathcal B}^{(n)}_{l_1} {\mathcal B}^{(n)}_{l_2} \cdot \frac{(-1)^{l_1+l_2}}{l_1! l_2!} b^{l_1} x^{l_2} \cdot e^{\imath \frac{\pi}{2} (l_2-l_1)} \end{eqnarray}

where

${\mathcal B}^{(n)}_{l} := \sum\limits_{j_1=l}^n (2(n-j_1)+1)^{(j_1-1)} j_1 (2(n-j_1)-1)!!$.

In the forth line we simplified the result and defined ${\mathfrak A}(x,b) := (x^{1-\theta} - b^{1-\theta})/(1-\theta)$ and in the fifth line we completed the argument of the exponential to a square. In the sixth line we decomposed a part of the integrand into partial fractions and integrated term by term. Here by$\left( \zeta_\xi^{(n)}\right)_{\xi=1}^{2 n}$ we denote the roots of a polynomial ${\mathbb C} \ni z \rightarrow \sum\limits_{j=0}^n (n-j+1)^{(2 j)}(2j-1)!!/(2^j j!)z^{2(n-j)} \in {\mathbb C} $. We have checked that for every $n\in {\mathbb N}$ those roots are not degenerate, i.e. all of them have a multiplicity one. We plot those roots below where the different colors correspond to different values of $n=0,\cdots,20$ (from violet to red respectively). We have:

b =.; num = 20;
(*Roots of the denominator.*)
cols = Table[ColorData["Rainbow", i/(num - 1)], {i, 0, num - 1}];
dens = Table[ 
   Sum[ Pochhammer[n - j + 1, 2 j] (b)^(2 (n - j))  (2 j - 1)!!/(
     2^j j!), {j, 0, n}], {n, 0, num}];
ListPlot[Transpose[{Re[#], Im[#]}] & /@ ((b /. NSolve[# == 0, b]) & /@
     Drop[dens, 1]), PlotStyle -> cols, AxesLabel :> {"Re", "Im"}, 
 PlotLabel :> "Zeros of BesselJ[n+1/2,z]^2+BesselY[n+1/2,z]^2"]

The final, seventh line, we used the identity ${\mathfrak C}^{(n,j)}_\xi + (-1)^j {\mathfrak C}^{(n,j)}_{2 n+1-\xi} = 0$ for $\xi=1,\cdots,n$ and $j=0,\cdots,2n+1$. Then we also used the fact that $\zeta^{(n)}_{2n+1-\xi} = \bar{\zeta^{(n)}_\xi}$, i.e. the $\xi$th and the $(2n+1-\xi)$th roots are complex conjugates.

Now the only thing that remains is to evaluate the improper integral, i,e, the term in the last line that is underlined. We refrain from writing it down explicitly in order not to obscure the whole picture but it is clear that that integral always exists; this follows form the fact that the roots are never purely real and negative; and the integral in question is expressed in terms of the error function the the hyperbolic cosine and sine integrals. This completes the derivation.

Below is a snippet of Mathematica code that verifies all the steps above, starting from the first line all the way down to the fifth. Here we go:

In[3523]:= (*The integrand*)
x =.; b =.; nu =.; Clear[f1]; Clear[f2]; Clear[FF];
f1[\[Nu]_, z_] := 
  Sum[((-1)^
       j (2 j + Abs[\[Nu]] + 1/2)!)/((2 j + 1)! (-2 j + Abs[\[Nu]] - 
         3/2)! (2 z)^(2 j + 1)), {j, 0, 
    Floor[(1/4) (2 Abs[\[Nu]] - 3)]}];
f2[\[Nu]_, z_] := 
  Sum[((-1)^
       j (2 j + Abs[\[Nu]] - 1/2)!)/((2 j + 0)! (-2 j + Abs[\[Nu]] - 
         1/2)! (2 z)^(2 j + 0)), {j, 0, 
    Floor[(1/4) (2 Abs[\[Nu]] - 1)]}];
(*For nu = 1/2 + n, with n being an integer, the function FF[] below \
is a multinomial in x and b of orders n,n respectively..*)
FF[nu_, x_, 
   b_] := (b x)^
    n ((-f2[nu, x] f1[nu, b] + f1[nu, x] f2[nu, b]) - 
     I (+f1[nu, b] f1[nu, x] + f2[nu, b] f2[nu, x]));
num = 4;
myasmpts := 0 < b < x && Element[n, Integers] && n >= 0;
mint0 = Simplify[
   ComplexExpand[
    Table[With[{nu = 1/2 + n}, 
      Im[BesselK[nu, -I x]/BesselK[nu, -I b]]], {n, 0, num}]], 
   Assumptions -> myasmpts];
mint = Table[
   With[{nu = 1/2 + n}, (
    BesselJ[nu, b] BesselY[nu, x] - BesselJ[nu, x] BesselY[nu, b])/(
    BesselJ[nu, b]^2 + BesselY[nu, b]^2)], {n, 0, num}];
mint1 = Table[With[{nu = 1/2 + n},
    (b/x)^(n + 1/2) (
         (Exp[+I (x - b)]) FF[nu, x, b]
       )/Sum[ 
       Pochhammer[n - j + 1, 2 j] (b)^(2 (n - j))  (2 j - 1)!!/(
        2^j j!), {j, 0, n}] ], {n, 0, num}];
mint1 = Simplify[ComplexExpand[Re[mint1]], Assumptions -> myasmpts];
(mint0 - mint) // PowerExpand // Simplify
(mint - mint1) // PowerExpand // Simplify
Out[3533]= {0, 0, 0, 0, 0}
Out[3534]= {0, 0, 0, 0, 0}

Answer 3

Here we are going to explain how do we invert the Laplace transform given in the body of the question, equation $(1)$. The idea is to expand the Laplace transform into powers of $(1-\theta)$ and then to invert term by term. The only problem is that even taking the limit $\theta \rightarrow 1_-$ is already hard because, in that case, both the index of the modified Bessel functions and their arguments go to plus infinity. Fortunately there is an asymptotic expansion that holds true in such a case, as given here, and we are going to apply it. We define $\nu: = (\mu_0 - 1/2)/(1-\theta)$ and $\delta(\alpha) := (2 \mu_0-1)^2 + 8 \alpha$, then $\eta_z := (1+z^2)^{1/2} + \log(z/(1 + (1+z^2)^{1/2}))$ and finally $U_k(p)$ as polynomials of degree $3 k$ in $p$ -- see the website in question for their precise definition--and here we go:

\begin{eqnarray} &&\frac{K_\nu\left[\frac{x^{1-\theta}}{1-\theta} \sqrt{2 \alpha} \right]} {K_\nu\left[\frac{b^{1-\theta}}{1-\theta} \sqrt{2 \alpha} \right]} = \\ && e^{-\nu \left[ \eta_{\frac{\sqrt{2 \alpha}}{\mu_0-1/2} \cdot x^{1-\theta}} - \eta_{\frac{\sqrt{2 \alpha}}{\mu_0-1/2} \cdot b^{1-\theta}} \right]} \cdot \frac{\left( 1 + \frac{2 \alpha}{(\mu_0-1/2)^2} \cdot b^{2(1-\theta)}\right)^{1/4}} {\left( 1 + \frac{2 \alpha}{(\mu_0-1/2)^2} \cdot x^{2(1-\theta)}\right)^{1/4}} \cdot \frac{ \left.\sum\limits_{k=0}^\infty (-1)^k \frac{U_k(p)}{\nu^k} \right|_{p = (1+ \frac{2 \alpha}{(\mu_0-1/2)^2} \cdot x^{2(1-\theta)})^{-1}} } { \left.\sum\limits_{k=0}^\infty (-1)^k \frac{U_k(p)}{\nu^k} \right|_{p = (1+ \frac{2 \alpha}{(\mu_0-1/2)^2} \cdot b^{2(1-\theta)})^{-1}} } = \\ && e^{-\frac{1}{2} \log(\frac{x}{b}) \sqrt{ \delta(\alpha) } } \sum\limits_{n=0}^\infty \left( 1-\theta \right)^n \sum\limits_{j=-(3n-1) 1_{n\ge 1} }^n [ \delta(\alpha) ]^{\frac{j}{2}} \cdot {\mathfrak A}^{(\mu_0,x,b)}_{n,j} \end{eqnarray} Here $\left( {\mathfrak A}^{(\mu_0,x,b)}_{n,j} \right)_{n=0,j=-(3n-1) 1_{n\ge 1} }^{\infty,n}$ are certain coefficients that depend on $(\mu_0,x,b)$ only but do not depend on $\theta$.

As always, I verify all the steps using a snippet of Mathematica code. This code constructs the quantity given in the second line from the top and then expands it in a series in powers of $(1-\theta)$. Having done that we check the normalization of all the modes and we verify that the $n=0$-mode integrates to unity and all the $n\ge 1$-modes integrate to zero!. Finally, for each value of $n=0,1,2$, we plot the coefficients in question along with the norms, i.e. the integrals over the half-axis, of each $n,j$-mode. Here we go:

(*Here, we invert the Laplace transform of the first hitting time \
distribution. In all what follows we use the asymptotic expansions \
from https://dlmf.nist.gov/10.41#ii.*)
p =.; z =.; nu =.; mu0 =.; th =.; eps =.; dd =.; a = (dd - (1 - 
      2 mu0)^2)/8;
myasms := th < 1 && mu0 > 1/2 && x > 0 && b > 0;
Us = {1, 1/24 (3 p - 5 p^3), 1/1152 (81 p^2 - 462 p^4 + 385 p^6), 
   1/414720 (30375 p^3 - 369603 p^5 + 765765 p^7 - 425425 p^9)};
Clear[eta]; eta[z_] := Sqrt[1 + z^2] + Log[z/(1 + Sqrt[1 + z^2])];
T3n = Sum[(-1)^k Us[[1 + k]]/nu^k, {k, 0, Length[Us] - 1}] /. 
     p :> 1/Sqrt[1 + z^2] /. z :> Sqrt[2 a]/(mu0 - 1/2) x^(1 - th) /. 
   nu :> (mu0 - 1/2)/(1 - th);
T3d = Sum[(-1)^k Us[[1 + k]]/nu^k, {k, 0, Length[Us] - 1}] /. 
     p :> 1/Sqrt[1 + z^2] /. z :> Sqrt[2 a]/(mu0 - 1/2) b^(1 - th) /. 
   nu :> (mu0 - 1/2)/(1 - th);
M = 3;
(T1=Exp[-nu (eta[Sqrt[2 a]/(mu0-1/2) x^(1-th)]-eta[Sqrt[2 

a]/(mu0-1/2) b^(1-th)])];)
res = FullSimplify[
   CoefficientList[
    Normal[Series[
       eta[Sqrt[2 a]/(mu0 - 1/2) x^(1 - th)] - 
        eta[Sqrt[2 a]/(mu0 - 1/2) b^(1 - th)], {th, 1, M + 1}, 
       Assumptions -> myasms]] /. th :> 1 - eps, eps]];
T1 = Exp[Collect[-(mu0 - 1/2)/(eps) (res.(eps)^Range[0, M + 1]), eps, 
    Simplify]];
(T2 =(( 1 + 2 a/(mu0-1/2)^2 b^(2-2 th))/(1 + 2 a/(mu0-1/2)^2 x^(2-2 

th)))^(1/4);)
T2 = FullSimplify[
   CoefficientList[
    Normal[Series[(( 
           1 + 2 a/(mu0 - 1/2)^2 b^(2 - 2 th))/(1 + 
            2 a/(mu0 - 1/2)^2 x^(2 - 2 th)))^(1/4), {th, 1, M}, 
       Assumptions -> myasms]] /. th :> 1 - eps, eps]];
T2 = T2.(eps)^Range[0, M];
T3 = FullSimplify[
   CoefficientList[
    Normal[Series[T3n/T3d, {th, 1, M}, Assumptions -> myasms]] /. 
     th :> 1 - eps, eps]];
T3 = T3.(eps)^Range[0, M];
(Expand the Laplace transform in a power series in (1-th).)
res = Together[
   CoefficientList[
    Normal[Series[T1 T2 T3, {eps, 0, M}, Assumptions -> myasms]], 
    eps]];
a =.;
(Check normalization. The "ground state" integrates to unity and all 

the "excited states" integrate to zero!)
Table[Simplify[((x/b)^(1/2 - mu0) res[[1 + n]] /. 
     dd :> (2 mu0 - 1)^2 + 8 a /. a :> 0), 
  Assumptions -> x > 0 && b > 0 && mu0 < 1/2], {n, 0, Length[res] - 1}]
(Plot the coefficients A_{n,j})
mcfs = FullSimplify[
   Table[{Coefficient[res[[1 + n]]/E^(1/2 Sqrt[dd] (Log[b] - Log[x])),
       dd, j/2], (2 mu0 - 1)^(j)}, {n, 0, 
     Min[Length[res] - 1, 2]}, {j, If[n == 0, 0, -(3 n - 1)], n}]];
MatrixForm[#] & /@ mcfs

Now we invert the Laplace transform term by term and we give a closed form expression for the probability density function of the first hitting time $\tau_b := inf( t\ge 0| X_t < b)$. Here we go:

\begin{eqnarray} \rho_{\tau_b}(t) &=& \left( \frac{x}{b}\right)^{1/2-\mu_0} \sum\limits_{n=0}^\infty \left(1-\theta\right)^n \sum\limits_{j=-(3n-1)1_{n\ge 1}}^n {\mathfrak A}^{(\mu_0,x,b)}_{n,j} \cdot \underline{ \frac{1}{8} e^{-1/2 (\mu_0-1/2)^2 t} \cdot {\mathcal L}_s^{-1} \left[ e^{-1/2 \log(\frac{x}{b}) \sqrt{s}} \cdot s^{j/2}\right](\frac{t}{8})} \\ &=& \left( \frac{x}{b}\right)^{1/2-\mu_0} \cdot \sum\limits_{n=0}^\infty \left(1-\theta\right)^n \sum\limits_{j=-(3n-1)1_{n\ge 1}}^n {\mathfrak A}^{(\mu_0,x,b)}_{n,j} \cdot % \frac{1}{8} e^{-1/2 (\mu_0-1/2)^2 t} % \cdot \left( \right. \\ && \left. \frac{1}{\Gamma(-j/2) (t/8)^{1+j/2}} F_{1,1} \left[ 1 + \frac{j}{2}, \frac{1}{2}; - \frac{[\log(x/b)]^2}{2 t} \right] \right. \\ && \left. - \frac{\log(x/b)}{2 (t/8)^{3/2+j/2} \Gamma(-1/2-j/2)} \cdot F_{1,1} \left[ \frac{3}{2} + \frac{j}{2}, \frac{3}{2}; - \frac{[\log(x/b)]^2}{2 t} \right] \right. \\ && \left. \right) \\ % &=& \left( \frac{x}{b}\right)^{1/2-\mu_0} \cdot \sum\limits_{n=0}^\infty \left(1-\theta\right)^n \sum\limits_{j=-(3n-1)1_{n\ge 1}}^n {\mathfrak A}^{(\mu_0,x,b)}_{n,j} \cdot % \frac{1}{8} e^{-1/2 (\mu_0-1/2)^2 t} % \cdot \left( \right. \\ && % (\frac{t}{8})^{-j/2-1} 2^{-j-1} \frac{1}{\sqrt{\pi}} \cdot \exp\left( - \frac{[\log(\frac{x}{b})]^2}{2 t}\right) \cdot H_{j+1}\left[ \frac{\log(\frac{x}{b})}{\sqrt{2 t}}\right] \\ % && \left. \right) % % % \end{eqnarray} Note, in the last step above, we have used the this identity. Here $H_j(\cdot)$ are the Hermite polynomials.

Now we have taken $\theta = 0.9$ and then we have chosen randomly $x,b,\mu_0$ subject to $x > b >0$ and $\mu_0 < 1/2$ and we produced two plots. On the left we plotted the cumulative sum of the first three modes ($n=0,\cdots, 2$) in equation $(2)$ above (violet to red for ascending values of $n$) along with the Bromwich integral itself (black). On the right we showed the absolute value of the difference between equation $(2)$ being truncated at $n=2$ and the Bromwich integral itself.

(*Compute the pdf of the first hitting time as a power series in \
(1-th).Run the block (*Here, we invert the Laplace transform *) to \
generate the quantities res[[]] first.*)
x =.; b =.; mu0 =.; th =.; t =.;
mypdf = (x/b)^(1/2 - mu0) Table[
    Sum[mcfs[[1 + n]][[j + If[n == 0, 1, 3 n], 1]] 1/
        8 Exp[-1/2 (mu0 - 1/2)^2 t] (t/8) ^(-1 - 
        j/2) (Hypergeometric1F1[1 + j/2, 1/2, -(Log[x/b]^2/(2 t))]/
         Gamma[-(j/2)] - (
         Log[x/b]/2 Hypergeometric1F1[(3 + j)/2, 3/
           2, -( (Log[x/b]^2)/(2 t))])/((t/8)^(1/2)
           Gamma[-(1/2) - j/2])), {j, If[n >= 1, -(3 n - 1), 0], 
       n}] (1 - th)^n, {n, 0, Length[mcfs] - 1}];
(Do the plotting. Run the block (Compute the pdf of the first h) 

to generate mypdf[[]] first.)
g = 0;
{x, b} = RandomReal[{0, 2}, 2, WorkingPrecision -> 50]; If[b > x, 
 tmp = x; x = b; b = tmp;];
mu0 = RandomReal[{-2, 0}, WorkingPrecision -> 50];
th = 9/10;
nu = (1 - 2 mu0)/(2 (1 - th));
t = Array[10^# &, 50, {-2, 1}];
cols = Table[
  ColorData["Rainbow", i/(Length[mypdf] - 1)], {i, 0, 
   Length[mypdf] - 1}]; cols = Join[cols, {Black}];
vals = Accumulate[Table[mypdf[[1 + n]], {n, 0, Length[mypdf] - 1}]];
vals0 = NIntegrate[
      Exp[(I a + g) #] (x/b)^(1/2 - mu0) BesselK[nu, 
         x^(1 - th)/(1 - th) Sqrt[2 (I a + g)]]/
        BesselK[nu, 
         b^(1 - th)/(1 - th) Sqrt[2 (I a + g)]], {a, -Infinity, 
       Infinity}]/(2 Pi) & /@ t;
ids = DeleteCases[Flatten[Position[vals0, _?(! NumberQ[#] &), 1]], 0];
 vals0[[ids]] = NaN;
vals = Join[vals, {vals0}];
pl1 = ListLogLinearPlot[Transpose[{t, #}] & /@ vals, PlotRange :> All,
    PlotStyle :> cols, Joined :> True, PlotMarkers :> Automatic, 
   AxesLabel -> {"t", "pdf[t]"}];
pl2 = ListLogLogPlot[Transpose[{t, Abs[vals[[3]] - vals0]}], 
   Joined :> True, PlotMarkers :> Automatic, 
   AxesLabel -> {"t", "|pdf[t]-pdf0[t]|"}];
GraphicsGrid[{{pl1, pl2}}, 
 PlotLabel :> "x,b,mu0,th=" <> ToString[N[{x, b, mu0, th}]]]
(t=.;Table[NIntegrate[mypdf[[1+n]],{t,0,10}] + 

NIntegrate[mypdf[[1+n]],{t,10,Infinity}],{n,0,Length[mypdf]-1}])

As we can see for $\theta = 0.9$ truncating the expansion $(2)$ at the first three terms, $n \le 2$, is completely sufficient for practical purposes.