Let $\{X_n:n=0,1,2,\ldots\}$ be a Markov chain on state space $S=\{0,1,\ldots,N\}$ with transition probabilities $P_{ij}$, $i,j\in S$. Define for $i\in S$ the hitting time $$ \tau_i = \inf\{n>0: X_n = i\mid X_0=0\}, $$ where $\tau_i:=+\infty$ when $$\mathbb P\left(\bigcup_{n=1}^\infty \{X_n=i\mid X_0=0\}\right)=0.$$ How do we compute the distribution of $\tau_i$? I have found many results for computing $\mathbb E[\tau_i]$, but that is not what I want. I am looking for a closed form for $$ \mathbb P(\tau_i = n),\ n=1,2,\ldots. $$ I have surveyed the literature and mostly found results for processes with additional assumptions, such as birth-death, skip-free, stationary, ergodic, aperiodic, irreducible, reversible, etc. Is there a general result for the distribution of the hitting times without these assumptions?
If not, what about in the case where the transition matrix is upper-triangular, that is, $P_{ij}>0$ if and only if $j\geqslant i$? Such a chain is necessarily absorbing, and it is well-known how to compute the expected time to absorption, but I cannot find any results on the distribution. Any references would be greatly appreciated.
I'm not sure if the following could be called a "closed form" or not, but I doubt that there's any other expression that could be better described as such. If $\ P_{(i)}\ $ is the matrix obtained from $\ P\ $ by replacing its $\ i^\text{th}\ $ column with a column of zeroes, then $$ \mathbb{P}\left(\tau_i=n\left|X_0=0\right.\right) = \mathbf{e}_0^TP_{(i)}^{n-1}P \mathbf{e}_i\ . $$
Example:
If $\ P= \pmatrix{\frac{1}{2}& \frac{3}{10}&\frac{1}{5}\\ 0&\frac{3}{5}&\frac{2}{5}\\ 0&0&1}\ $, as the OP has asked to be treated in the comments, then \begin{align} P_{(0)}&=\pmatrix{0& \frac{3}{10}&\frac{1}{5}\\ 0&\frac{3}{5}&\frac{2}{5}\\ 0&0&1}\ ,\\ P_{(1)}&=\pmatrix{\frac{1}{2}&0&\frac{1}{5}\\ 0&0&\frac{2}{5}\\ 0&0&1}\text{ , and }\\ P_{(2)}&=\pmatrix{\frac{1}{2}& \frac{3}{10}&0\\ 0& \frac{3}{5}&0\\ 0&0&0}\ . \end{align} For small examples like this, it is possible to obtain expressions for $\ \mathbb{P}\left(\tau_i=n\left|X_0=0\right.\right)\ $ that are unequivocally describable as "closed form".
For $\ n=1,2,\dots\ $ we have \begin{align} P_{(0)} ^n &=\pmatrix{0&\frac{1}{2} \left(\frac{3}{5}\right)^n& \frac{1}{2} \left(1-\left(\frac{3}{5}\right)^n\right)\\ 0&\left(\frac{3}{5}\right)^n& 1-\left(\frac{3}{5}\right)^n\\ 0&0&1}\ ,\\ P_{(1)}^n&=\pmatrix{\frac{1}{2^n}&0& \frac{2}{5}\left(1-\frac{1}{2^n}\right)\\ 0&0&\frac{2}{5}\\ 0&0&1}\\ P_{(2)}^n&=\pmatrix{\left(\frac{1}{2}\right)^n& 3 \left(\left(\frac{3}{5} \right)^n-\left(\frac{1}{2} \right)^n\right)&0\\ 0& \left(\frac{3}{5} \right)^n&0\\ 0&0&0}\ . \end{align} Substituting these into the above formula for $\ \mathbb{P}\left(\tau_i=n\left|X_0=0\right.\right)\ $, we get \begin{align} \mathbb{P}\left(\tau_0=n\left|X_0=0\right.\right)&=\cases{\frac{1}{2}&for $\ n=1\ $\\ 0 &for $\ n=2,3,\dots\ $, or\\ \frac{1}{2}&for $\ n=\infty\ $.}\\ \mathbb{P}\left(\tau_1=n\left|X_0=0\right.\right)&=\cases{\frac{3}{10}\left(\frac{1}{2}\right)^{n-1}&for $\ n=1,2,\dots\ $, or\\ \frac{2}{5}&for $\ n=\infty\ $.}\\ \mathbb{P}\left(\tau_2=n\left|X_0=0\right.\right)&=2\left(\frac{3}{5} \right)^n-\left(\frac{1}{2} \right)^{n-1}\ . \end{align} Note that the condition given by the OP for $\ \tau_i=\infty\ $ is incorrect. The event \begin{align} \left\{\tau_i=\infty\right\}&=\Omega\setminus\bigcup_{n=1}^\infty\left\{\tau_i=n\right\}\\ &= \Omega\setminus\bigcup_{m=1}^\infty\left\{X_m=i\right\} \end{align} has probability $\ 1-\mathbb{P}\left(\bigcup_{m=1}^\infty\left\{X_m=i\right\}\left|\,X_0=0\right.\right)\ $. Thus, $\ \mathbb{P}\left(\bigcup_{m=1}^\infty\left\{X_m=i\right\}\left|\,X_0=0\right.\right)=0\ $ means that $\ \tau_i=\infty\ $ with probability $1$. That is not the case for any of the $\ \tau_i\ $ in this example. Though $\ \tau_0\ $ and $\ \tau_1\ $ are infinite with positive probability, that probability is not $1$ in either case.