Given $\lbrace Y_i\rbrace$ a non-homogenous Poisson process with mean density $\theta y^{-1}e^{-y}$ where $y>0$ $(\theta>0)$. I.e., the number of points of $\lbrace Y_i\rbrace$ in $(a,b)$ with $0<a<b\le \infty$ has a Poisson distribution with mean $\int_a ^b \theta y^{-1}e^{-y}dy$.

How can the density of $Y_j$ be expressed as $$\theta y^{-1}\frac{(\theta E_1(y))^{j-1}}{(j-1)!}e^{-\theta E_1(y)}, y>0 \tag{$*$}$$

where $E_1(y)$ is the exponential integral $$E_1(y)=\int_y ^\infty \frac{e^{-w}}{w}dw\,\,?$$ I found this in page 3 of http://www.stats.ox.ac.uk/~griff/pd.pdf

I know how to find the number of $Y_j$ in a given interval: $$\mathbb{P}(N(b)-N(a)=n)=\frac{(\int_a ^b \theta y^{-1}e^{-y}dy)^n \exp(-\int_a ^b\theta y^{-1}e^{-y}dy )}{n!}$$

Although this formula has some similarities to the desired formula, I am not sure how to derive $(*)$.