I am trying to understand the basics of continuous time random walks, and this formula has no explanation as far as I can find:
$$\hat{p}_{n}(s) = \hat{\rho}(s)^n\cdot\frac{1-\hat{\rho}(s)}{s}$$
Where $p_{\hphantom{.}n}(t)$ is the probability of $n$ jumps in $t$ time (meaning that the $n$th jump has occured before time $t$ and the $n+1$-th will occur after time $t$), $\rho(t)$ is the PDF of the wait time distribution between jumps, and the hat above them denotes each one's respective Laplace transform.
I have seen an article that managed to develop $\hat{p}_{n}(s)$ into
$$\mathbb{E}\bigg[e^{-s t_n}\cdot\frac{1-e^{-s\tau_{n+1}}}{s}\bigg]$$
Where $t_n$ is the time of the $n$th jump and $\tau_{n+1}$ is $t_{n+1} - t_n$. But then it immediately said that's equal to the aforementioned formula with no explanation, and I couldn't find anything that could explain the equality earlier in the article, nor any citations.
If it's of any help, I also found another article that wrote it was "because of the convolution property of Laplace transforms".
Can someone explain why this is?
Not sure if this qualifies as an answer since I'm self-taught when it comes to Laplace transforms - but I think I've come across something similar whilst trying to handle the following adjacent problem:
Consider a Poisson-distributed variable $B(t)$ denoting a random event, so that the probability of $n$ event arrivals in any interval of length $t$ is $\mathrm{Pr}\left[B(t)=n \right]=\frac{(\lambda t)^n}{n!}e^{-\lambda t}, \ n=0,1,2,\dots$. Taking the Laplace transform of this expression you get the following: $$ \begin{align} \mathcal{L}\left\{ \mathrm{Pr}\left[B(t)= n \right] \right\} &= \mathcal{L}\left\{\frac{(\lambda t)^{n}}{n!}e^{-\lambda t} \right\}\\ &=\int_{0}^{\infty} e^{st}\frac{(\lambda t)^{n}e^{-\lambda t}}{n!} dt \\ &=\frac{(\lambda t)^{n}}{n!}\int_{0}^{\infty} t^{n}e^{-t(s+\lambda)}dt \\ &=\frac{(\lambda t)^{n}}{n!}\left[ \frac{n!}{(\lambda+s)^n} \frac{1}{(\lambda+s)} \right] \\ &=\left(\frac{\lambda}{\lambda+s} \right)^{n+1}\frac{1}{\lambda} \\ &=\left(\frac{\lambda}{\lambda+s} \right)^{n}\frac{1}{\lambda+s} \end{align} $$
Now, you might already know that the Poisson and negative exponential are related Wagner (1975) p.860. So much so that the previous expression can be written in terms of the transform of the negative exponential density function - what you call $\rho(t)$ in your query i.e. $\rho(t) = \lambda e^{-\lambda t}$. This becomes evident once you take the Laplace transform of the negative exponential density function:
$$ \begin{align} \mathcal{L}\left\{ \lambda e^{-\lambda t} \right\} &= \int_{0}^{\infty}e^{-st}\lambda e^{-\lambda t}dt \\ &= \lambda \int_{0}^{\infty}e^{-t(s+\lambda)} \\ &= \lambda \left[ \frac{1}{s+\lambda} \right] \\ &= \frac{\lambda}{s + \lambda} \\ &= \widehat{\rho}(s) \end{align} $$
Also notice that, based on the above: $ 1-\mathcal{L}\left\{ \lambda e^{-\lambda t} \right\} = \frac{s}{\lambda + s} $
This means that we can rewrite the initial expression in such a way that it yields the right-hand side of the equation you're inquiring about i.e. $$ \begin{align} \mathcal{L}\left\{ \mathrm{Pr}\left[B(t) = n \right] \right\} &= \mathcal{L}\left\{ \lambda e^{-\lambda t} \right\}^{n} \frac{1}{s}\left( 1 - \mathcal{L}\left\{ \lambda e^{-\lambda t} \right\} \right) \\ &= \widehat{\rho}(s)^n \frac{1}{s}\left( 1 - \widehat{\rho}(s) \right) \end{align} $$
Finally, regarding your point on the convolution properties of the Laplace transform: Muth (1977) p.196-197 shows that the density function of the sum of two independent random variables equals the convolution of the density functions of the two respective random variables.
For example, if $X$ and $Y$ are two independent random variables with probability density function $f_{X}(t)$ and $f_{Y}(t)$, respectively, then their sum $Z = X + Y$ has probability density function $f_{Z}(t)$ whose Laplace transform is $\mathcal{L}\left\{ f_{Z}(t)\right\}=E\left[ e^{-s(X+Y)}\right]=E\left[ e^{-sX}\right]E\left[e^{-sY}\right]= \mathcal{L}\left\{ f_X(t) \right\} \mathcal{L}\left\{ f_Y(t) \right\}=\tilde{f}_{X}\tilde{f}_{Y}$. Therefore $ f_{X+Y}(t) = \mathcal{L}^{-1}\left\{ \tilde{f}_{X}(s)\tilde{f}_{Y}(s) \right\}= f_X(t) \ast f_Y(t)$. This result builds on the fact that the expected value of a function of a random variable X is $E[g(X)]=\int_{0}^{\infty}g(t)f_X(t)dt$, and that for $g(X)=e^{-sX}$ we can write $E[e^{-sX}]=\mathcal{L}\left\{f_X(t) \right\}$.
In the specific case in which $X$ and $Y$ follow the negative exponential densifty function i.e. $f_X(t)=f_Y(t)=\lambda e^{-\lambda t}$, you have that the Laplace transform of their sum $Z = X + Y$ is as follows: $$ \begin{align} \mathcal{L}\left\{f_Z(t) \right\} &= E\left[ e^{-sX}\right]E\left[e^{-sY}\right] \\ &=\left( \frac{\lambda}{s + \lambda} \right)^2 \end{align} $$ which is just $\widehat{\rho}(s)^n$ when $\rho(t) = \lambda e^{-\lambda t}$ and $n=2$
I hope the above is somewhat useful