I have read the paper Combinatorics of geometrically distributed random variables: Left-to-right maxima. In the paper, the largest order statistic $X_{n:n}$ (i.e., $\max\{X_1,X_2,\ldots,X_n\}$) is analyzed, where $X_i$ are iid geometrically distributed rvs and $\Pr \{X_i=x \} = pq^{x-1}$. $\mathbb E[X_{n:n}]$ is analyzed by the use of Rice's method.
I am interested in $\mathbb E[X_{k:n}]$, i.e., the expected value of the $k$th order statistic. I googled and found only the paper Combinatorics of Geometrically Distributed Random Variables: Value and Position of the $r$th Left–to–right Maximum. But the $r$th left–to–right maximum is different from the $r$th order statistic.
Moreover, letting $\Delta = n-k$, I derived an formula for $\mathbb E[X_{k:n}]$ as \begin{align} & \mathbb E[X_{k:n}] = 1+\sum_{m=0}^{k-1} {n \choose m} \sum_{r=0}^m {m \choose r} (-1)^r \cdot \frac{q^{n+r-m}}{1-q^{n+r-m}} \\ = & 1+\sum_{m=n-k+1}^n (-1)^{1+m-(n-k)} \cdot \frac{m-(n-k)}{m} {m \choose m-(n-k)} {n \choose m} \cdot \frac{q^m}{1-q^m} \\ = & 1+\sum_{m=1}^n (-1)^{1+m-\Delta} {n \choose m} {m-1 \choose \Delta} \cdot \frac{q^m}{1-q^m}\\ =& 1+\sum_{m=1}^n {n \choose m} (-1)^{m+1} {\Delta -m \choose \Delta} \cdot \frac{q^m}{1-q^m} . \tag{1} \end{align}
Substituting $k=n$ into (1), We can see that \begin{align} \mathbb E[X_{n:n}] = 1+ \sum_{m=1}^n (-1)^{m+1} {n \choose m} \cdot \frac{q^m}{1-q^m}. \tag{2} \end{align}
Compared to (2), which is directly applicable to Rice's method, the term in (1) ${\Delta -m \choose \Delta} \cdot \frac{q^m}{1-q^m}$ is not of polynomial growth. Thus I don't know how to apply Rice's method to do the asymptotic analysis.
Use $Geom(p=1-q)\overset{d}=\lfloor\mathcal{Exp}(\lambda=-\log q)\rfloor$ (set $\lambda=1$). You immediately get $$-1<E(X^g_{k:n})-E(X^e_{k:n})<0$$ which yields first-order asymptotics for $k=n-o(n)$ (e.g. $k=n$).
For $k=o(n)$ (e.g. $k$ fixed) we get: $$E(X^g_{k:n})= E(\lfloor X^e_{k:n}\rfloor|X^e_{k:n}\ge 1)P(X^e_{k:n}\ge 1)=(1+O(\frac 1 n)) P(X^e_{k:n}\ge 1)\tag 0$$ Let $Z=\sum_{i=1}^k \mathcal{Exp}(\lambda_i)$. It can be easily shown by induction and symmetry that $$P(Z\ge x)=\sum_{i=1}^k(e^{-\lambda_ix}\prod_{j\ne i}\frac {\lambda_j}{\lambda_j-\lambda_i})\tag 1$$
Plugging in $\lambda_i=n-i+1$ for $X^e_{k:n}$ the above expression simplifies to $$P(X^e_{k:n}\ge x)=e^{-(n-k+1)x}\binom{n}{k-1}(1-e^{-x})^{k-1}(1+O(\frac 1 n))\tag 2$$
This matches the fact that $X^e_{k:n}$ lies between $Erlang(k,n)$ and $Erlang(k,n-k+1)$ and hence $e^{-n}\frac {n^{k-1}}{(k-1)!}$ and $e^{-n+k-1}\frac {(n-k+1)^{k-1}}{(k-1)!}$ are lower and upper asymptotic boundaries of $P(X^e_{k:n}\ge 1)$.
ETA: From $$f_{X^e_{k:n}}(x) = n {n-1\choose k-1} (1-e^{-x})^{k-1} e^{-(n-k+1)x}$$ for $\limsup \frac{k-1}n<1-e^{-1}=F_{\lambda}(1)\approx 0.62$ Laplace methods yields: $$P(X^e_{k:n}\ge 1) \sim E(X^e_{k:n};X^e_{k:n}\ge 1)\sim \binom {n-1}{k-1} {e^{-(n-k+1)}(1-e^{-1})^{k-1}}\frac 1{1-\frac{k-1}{n(1-e^{-1})}}$$ and same is true for $\lfloor X^e_{k:n}\rfloor$ and yields first-order asymptotics for a much wider range of $k$ - including those proportional to $n$. More generally, for $F_{\lambda}^{-1}(\frac k n\to\alpha<1)=Q_\alpha$ not a positive integer, $$E(X^g_{k:n})-\lfloor Q_{\alpha}\rfloor=E(\lfloor X^e_{k:n}\rfloor)-\lfloor Q_{\alpha}\rfloor\sim P(X^e_{k:n}\ge \lfloor Q_{\alpha}\rfloor+1)-P(X^e_{k:n}< \lfloor Q_{\alpha}\rfloor)$$ and the last two integrals can be similarly derived. For $Q_\alpha$ a positive integer, the situation is a little bit more involved and depends on behaviour of $n^{1/2}(\frac {k-1}n-\alpha)$ since $n^{1/2}(X^e_{k:n}-Q_\alpha)\to N(0,\frac {\alpha}{1-\alpha})$.