Probability of rolling a 1 before you roll two 2's, three 3's, etc

This is a question my father-in-law asked. I found it interesting but haven't been able to answer it.

Suppose you have an $N$ sided die. You conduct an experiment as follows: roll the die until the first time you have rolled one 1, two 2's, three 3's etc. For example, you roll 2,3,4,3,6,5,3 the experiment ends because you have rolled three 3's. Of course, the most times you can roll the die is $1+\sum_{k=1}^{N-1}k = 1 + \frac{(N-1)N}{2}$. This is because the maximum number of rolls you can do is to roll one 2, two 3's, three 4's and so on, and then one more roll to complete the experiment. In particular, the experiment always ends in a finite number of rolls.

Observe that the experiment ends when you roll a $1$. If you have a one-sided die the experiment always ends after the first roll and so you end on $1$ with probability 1. If you have a two sided die (i.e. a coin) then the probability of rolling a $1$ on the first roll is $1/2$; the probability of rolling a 2 first and a 1 second is $1/4$ and so the probability of ending the experiment on a $1$ is $3/4$.

Let $P_N$ be the probability that you end the experiment on a $1$. What is $\lim_{N\to\infty}P_N$?

$P_N$ is decreasing and of course $0\leq P_N$ so the limit exists. I tried taking the limit along ``easy'' sequences like $2^N$ hoping the the extra structure would lend itself to easier analysis, but I couldn't.

I imagine I am not the first to ask this question, but I couldn't find it on MSE any other place. Any references to similar questions are also appreciated.

I doubt that a "closed form" for their limit will be forthcoming, but the probabilities can be written $$P(n) = \sum_{k=0}^{T_{n-1}}\frac{a(k,n)}{n^{k+1}}$$ where $T_{n-1}=\frac{(n-1)n}{2}$ is the $(n-1)$th triangular number, and $a(k,n)$ is the number of $k$-long words that can be formed using $k$ elements from the string $S_n=1^02^13^2...n^{n-1}$ such that no element is used more often than it appears in $S_n$. (The string notation $c^i$ denotes $i$ concatenated copies of the symbol $c$.)

Now, $a(k,n)$ is readily computed, as it is just $k!$ times the coefficient of $x^k$ in the following exponential generating function:

$$\left(1+x\right)\,\left(1+x+\frac{x^2}{2!}\right)\,\left(1+x+\frac{x^2}{2!}+\frac{x^3}{3!}\right)\,...\,\left(1+x+\frac{x^2}{2!}+...+\frac{x^{n-1}}{(n-1)!}\right). $$

(A simple explanation of this type of generating function is in the accepted answer to the MSE question 6-letter permutations in MISSISSIPPI.)

Computations produce the following values:

n     P(n)
--    -----------------
 1    1.000000000000000
 2    0.750000000000000
 3    0.703703703703704
 4    0.692138671875000
 5    0.689010688000000
10    0.687835012279842
15    0.687833655173612
20    0.687833654277303
25    0.687833654276930
30    0.687833654276930

The online Inverse Symbolic Calculators found no closed form for this apparent limit of 0.687833654276930....

For reference, here's the (Sage) program I used:

def P(n):
    return sum( a(k,n)/n^(k+1) for k in [0..n*(n-1)/2] )

def a(k,n):
    R.<x> = PowerSeriesRing(ZZ)
    f = prod( sum( x^j / factorial(j) for j in [0..i] ) for i in [0..n-1] )
    return factorial(k)*f[k]

for n in [1..50]:
    print n, P(n).n()