Let the interarrival times be identically $1$ and the service rate $\mu>1$ (so that $1/\mu<1$). Let $B$ be the distribution of the number of service completions that occur between the arrival of two customers, then $$b_j:=\mathbb P(B=j) = \int_0^\infty e^{-\mu t}\frac{(\mu t)^j}{j!}\,\mathsf dA(t) = e^{-\mu}\frac{\mu^j}{j!}, $$ where $A(t)=\mathsf 1_{[1,\infty)}(t)$ is the distribution function of the interarrival time. Let $\{X_n:n\in\mathbb N_0\}$, where $X_n$ is the number of customers in the system immediately before the $n^{\mathrm{th}}$ arrival, then the transition probabilities are given by $$ P_{ij} = \begin{cases} b_{j-i+1},& j\geqslant i-1\\ 1-\sum_{k=0}^i b_k,& j=0. \end{cases} $$ From $\pi=\pi P$ we have $$\pi_k = \sum_{i=k-1}^\infty \pi_i P_{ik} = \sum_{i=k-1}^\infty \pi_i b_{i+1-k} =\sum_{i=k-1}^\infty \pi_i e^{-\mu} \frac{\mu^{i+1-k}}{(i+1-k)!}, $$ which is a countably infinite system of linear equations with no evident solution.
If we assume that $\pi_k = (1-\xi)\xi^k$ for some $\xi\in(0,1)$, then dividing by $(1-\xi)\xi^{k-1}$ yields $$\xi = \sum_{i=k-1}^\infty \xi^{i+1-k} e^{-\mu}\frac{\mu^{i+1-k}}{(i+1-k)!}=e^{-\mu}\sum_{j=0}^\infty \frac{(\mu\xi)^j}{j!} = \exp\left(-(\mu-\mu\xi)\right). $$ Now, this is precisely the equation to determine the extinction probability in a Galton-Watson process with a Poisson offspring distribution. This makes sense intuitively, as this is fundamentally the same process as the D/M/1 queue (conditioned on there being one customer in the system at time zero and killed when the system becomes empty). From here it follows from the Banach fixed point theorem that for an arbitrary $\xi_0\in(0,1)$, the sequence $\xi_{n+1} = \exp(-(\mu-\mu\xi_n))$ satisfies $\lim_{n\to\infty}\xi_n = \xi$, yielding an iterative scheme to compute $\xi$ numerically.
The only problem with this argument is that the assumption $\pi_k=(1-\xi)\xi^k$ is rather arbitrary. It seems like there is a more fundamental relationship between the G/M/1 queue and the Galton-Watson process which one could use to show that the stationary distribution is geometric in form. I am not sure how to make this rigorous though, and would appreciate suggestions or insight into this problem.