Consider a continuous time random walk which steps by an exponentially random length $l$ at time-points dictated by a Poisson process of rate $\omega$:
$$ x(t) = \sum_{k=1}^{N(t)} l_k.$$
The number of steps $N(t)$ is distributed as $[(\omega t)^n/n!]e^{-\omega t}$. The Master equation for the probability of being at position $x$ by time $t$ is
$$ \partial_t p(x,t) = -\omega p(x,t) + \omega \int_0^\infty dl p(x-l,t)h(l),$$ where $h(l)=\lambda e^{-\lambda l}$ is the step length distribution. How can I find this $p(x,t)$? I expect it's proportional to the modified Bessel function of the first kind $I_1(2\sqrt{\lambda x\omega t})$ from several papers, but these do not explain the derivation. A more general case is described in Cox and Miller (1965) (eq. 157), but they do not evaluate a closed form for $p(x,t)$.
So far, I have tried Fourier transform methods without much success. Fourier transforming both sides gives
$$\partial_t \phi(s,t) = -\omega \phi(s,t) + \omega \phi(s,t)\hat{h}(s),$$ where $\hat{h}(s) = \lambda/(\lambda-is)$ is the FT of $h(l)$ and $\phi(s,t)$ is the characteristic function of $p(x,t)$.
Solving this and taking the inverse transform (using initial conditions $p(x,0)=\delta(x)$ or $\phi(s,0)=1$) gives
$$ p(x,t) = \frac{1}{2\pi} \int_0^\infty ds \exp\Big\{-isx + \frac{i \omega ts}{\lambda - i s} \Big\}, $$ which I can't get much traction on.
Any suggestion to proceed? I tried substitutions to arrive at an integral representation of a bessel function, but no progress yet. I can arrive at forms in integral tables but without the proper limits. Any approaches are appreciated. Thanks!
A better way to do this problem than the one discussed earlier in the comments is to marginalize over $N(t)$, meaning that you write
$$p(x,t)=f_{X(t)}(x)=\sum_{n=0}^\infty f_{X(t)}(x \mid N(t)=n) P(N(t)=n).$$
The reason this method is good is that $f_{X(t)}(x \mid N(t)=n)=\lambda^n \frac{x^{n-1}}{(n-1)!} e^{-\lambda x}$ and $P(N(t)=n)=e^{-\omega t} \frac{(\omega t)^n}{n!}$ are both known functions. Here we must understand $\frac{x^{-1}}{(-1)!} e^{-\lambda x}$ as $\delta(x)$; this makes sense because the resulting term represents the probability that no jumps happen, i.e. the particle stays at $x=0$ all the way until time $t$.
Plugging that in gives
$$p(x,t)=e^{-\omega t} \delta(x) + e^{-\omega t} e^{-\lambda x} \lambda \omega t \sum_{n=0}^\infty (\lambda x \omega t)^n \frac{1}{n!} \frac{1}{(n+1)!}.$$
As you mentioned in the OP you can now write this in a way like you expected, because $I_1(2\sqrt{\lambda x \omega t})=\sqrt{\lambda x \omega t} \sum_{n=0}^\infty \frac{(\lambda x \omega t)^n}{n!(n+1)!}$ from https://dlmf.nist.gov/10.25#E2 Therefore
$$p(x,t)=e^{-\omega t} \delta(x) + e^{-\omega t} e^{-\lambda x} \sqrt{\frac{\lambda \omega t}{x}} I_1 \left ( 2 \sqrt{\lambda x \omega t} \right ).$$