Generations until fixation: A nontrivial generalization of a dice convergence problem In spite of its "recreational" aspect, this question appears to me to be research-level and (I hope) clearly formulated and tagged.

Edit 4/4/20: You can find a related question with the same motivator as posed by Greg Egan:


Background:
Here is a relatively straightforward question from a recent 538 Riddler:

Last week, you started with a fair 6-sided die and rolled it six times, recording the results of each roll. You then wrote these numbers on the six faces of another, unlabeled fair die. For example, if your six rolls had been 3, 5, 3, 6, 1 and 2, then your second die wouldn’t have had a 4 on it; instead, it would have two 3s.
Next, you rolled this second die six times. You took those six numbers and wrote them on the faces of yet another fair die, and you continued this process of generating a new die from the previous one.
Eventually, you’d have a die with the same number on all six faces. What was the average number of rolls it would take to reach this state?

The link above contains an exact solution, but there is also a follow-up with some work done, yet which lacks a closed form solution.
My Question:

What is the precise answer for the version of this question in which you start with a fair $N$-sided die?

The 538 piece above includes a link to work from Angela Zhou, who also tweets:

 A: A similar idea was considered by Motoo Kimura and Tomoko Ohta in The Average Number of Generations Until Fixation of a Mutant Gene. They used a diffusion model, which should correspond to looking at the behavior as $N \rightarrow \infty$. The relevant line should be:

NOW, in a  population consisting of N individuals, if we assume that each mutant gene is represented only once at the moment of its  occurrence, p = 1/ (2N), and from formula (14), the average number of generations until fixation of a neutral mutation becomes
  $\bar{t}_1(\frac{1}{2N}) = -8 N N_e (1 - \frac{1}{2N}) \log_e(1 - \frac{1}{2N})$

Our $N$ is Kimura and Ohta's $2N$; if I understand correctly, their $N_e$ should equal $N$ in our situation ($N_e$ denotes an "effective" population size where individuals act differently; here, that isn't true). Correspondingly, our limiting behavior should be:
$-2 N^2 (1 - \frac{1}{N}) \ln(1 - \frac{1}{N}) = 2 N (N - 1) \sum_{i = 1}^\infty \frac{1}{i N^i} = 2 \left(\sum_{i = 1}^\infty \frac{1}{i N^{i - 2}} - \sum_{i = 1}^\infty\frac{1}{i N^{i - 1}}\right)$
$= 2N + 2\sum_{i = 0}^\infty \frac{1}{N^i} (-\frac{1}{i + 1} + \frac{1}{i + 2}) = 2N - 2\sum_{i = 0}^\infty \frac{1}{(i + 1)(i + 2) N^i}$
$ = 2N - 1 - \frac{1}{3N} - O(\frac{1}{N^2})$
This is subject to the diffusion model being accurate to this degree of approximation, which I have no specific support for. I wouldn't be surprised to see the logarithmic term in Angela Zhou's answer show up when considering the effect of the first "round". 
Tracing back the references a bit, I found this paper, which may be more obviously relevant (as well as already being written for mathematicians directly).
A: I can propose a simple heuristic but it doesn't go as far as Angela Zhou conjecture.
For each face the probability that it appears $k$ times in $N$ toss is $\frac{N!}{k!(N-k)!}\frac{(N-1)^{n-k}}{N^N}$. For $N$ large this can be approximate by the Poisson distribution $\frac{1}{k!}e^{-1}$. 
As repeating the process, for any fixed number the number of face of that number should behave rougthly like a Galton Watson process as each face are almost independent. 
So my claim is that this problem should behave for $N$ large as $N$ iid Galton Watson process with poisson distribution:
$Z_{n+1}:=\sum_{i=1}^{Z_n}P_{i,n}$ with $P_{i,n}$ iid Poisson variable. Considering the generating function
$$\mathbb{E}(s^{Z_1})=\sum_{k}\frac{s^k}{k!}e^{-1}=e^{s-1}$$we get that $\mathbb{E}(s^{Z_n})=u_n(s)$ where $u_0(s)=s$ and for all $n$ $u_{n+1}(s)=e^{u_n(s)-1}$.
In particular as $\mathbb{P}(Z_n=0)=u_n(0)$, the proportion of numbers that have disappeared at time $n$ is about $Nu_n(0)$, which are for $n=1:$ $Ne^{-1},$ $n=2: Ne^{e^{-1}-1},$ $n=3: Ne^{e^{e^{-1}-1}-1},\cdots$). 
As $u_n(0)\rightarrow 1$, we also have that $$u_{n+1}(0)-1=u_n(0)-1+\frac{1}{2}(u_n(0)-1)^2+\mathcal{O}((u_n(0)-1)^3)$$ which allow one to get the behavious $(u_n(0)-1)\sim-\frac{2}{n}$ (This assuming $N(u_n(0)-1)=1$ gives the $2N$ expected rolls).
A: This is the simplest (discrete) case of Kingman's coalescence. See e.g. https://arxiv.org/abs/0809.4233  and the explanations and references there.
The relation to the process considered there can be seen as follows:
view the results of the $N$ throws of the first round as value table $(f_1(1),\ldots,f_1(N))$ of a random mapping $f_1$.
Imagine that you write the pairs $(i,f_1(i))$ on the faces of the unlabeled die (instead of just the values $f_1(i)$, as you do). 
Then the results of the $N$ throws of the second round may be viewed as producing the value table $\big((f_2(1),f_1(f_2(1)),\ldots,(f_2(N),f_1(f_2(N))\big)$,
(where $f_2$ is a random mapping, and independent of $f_1$), 
and so on, and you are interested in the "coalescence time" $T$ of this process, i.e. the  smallest $t$ for which the composition $f_1\circ\ldots\circ f_t$ of independent (uniform)
$N$ to $N$ mappings becomes constant. 
