I have a Markov Chain of $N$ states. Such states represent the energy levels in a molecule.

The states' connectivity is as follows:

States $j\in\{0,\ldots,N\}$ transition to $k\in\{\max(j-M,0),...,\min(j+M,N)\}\setminus j$, with $N\gg M$. The resulting transition matrix $\mathbf{P}$ has known terms $P_{jk}$, and it is a hollowed band matrix where detailed balance is followed.

States $l\in\{L,...,N\}$ ALSO transition outside the Markov Chain, with rates $\mu_l>0$. Every state $l$ retains transitions to the $k$ states described in the previous point. For instance, state $L+1$ has both $j$ and $l$ transitions.

I am interested in recovering the stationary distribution $\tilde\pi_j$ of such Markov Chain.

If $\mu_l=0$, the Markov Chain has stationary distribution $\pi_j$ known from detailed balance and probability closure condition: $$ \pi_kP_{kj}=\pi_{j}P_{jk},\quad \sum_j\pi_j=1 $$ writing $\pi_j=\pi_0P_{0j}/P_{j0}$ and substituting into $\sum_j\pi_j=1$, one obtains $\pi_0$ and all the remaining values.

Now, when $\mu_l>0$, the stationary distribution $\tilde\pi_j$ is different. To tackle the problem, I have assumed that the $l$ states go to an absorbing state $J=N+1$, thus the transition matrix gains an absorbing state (one row and one column) with: $$ P_{lJ}=\mu_l,\quad P_{(j\setminus l)J}=0,\quad P_{Jl}=0,\quad P_{J(j\setminus l)}=0,\quad P_{JJ}=1 $$ But from now, to the best of my knowledge, I do not know how to proceed. I have the feeling that by introducing $J$, I am rendering the Markov Chain not ergodic anymore (as all states are unreachable from $J$). I have the feeling that $\tilde \pi_j=\delta(j-J)$, having $\delta(x)$ as the Dirac Delta Function.

I have read various sources (Snell, Montroll, Othmer) on Markov Chains and queuing theory, but all examples were not as specific as the one I have described above.