[EDITED mostly to report on the answer by Kevin Costello (and to improve the gp code at the end)]
I thank Nicolas Dupont for the following question (and for permission to disseminate it further):
[Nicely answered by Kevin Costello: $M_N$ is asymptotic to $e \log N$ as $N \rightarrow \infty$. Moreover, the maximal multiplicity is within $o(\log N)$ of $e \log N$ with probability $\rightarrow 1$. I don't recall any other instance of a naturally-arising asymptotic growth rate of $e \log N$...]
[Looks like this was a red herring: "One might guess" plausibly that $M_N - H_N$ is dwarfed by $H_N$ for large $N$, but by Kevin Costello's answer $M_N - H_N$ asymptotically exceeds $H_N$ by a factor $e - 1$, and that factor is more complicated than $\lim_{N\rightarrow\infty} M_N/H_N = e$, so analyzing the difference $M_N-H_N$ is likely not a fruitful approach.]
@ For each $N>1$, the expected value of the maximum count is given by the convergent $(N-1)$-fold sum $$ M_N = \sum_{a_1,\ldots,a_{N-1} \geq 1} N^{-\!\sum_i a_i} (\sum_i a_i)! \frac{\max_i a_i}{\prod_i a_i!}. $$$$ M_N = \sum_{a_1,\ldots,a_{N-1} \geq 1} N^{-\!\sum_i \! a_i} \Bigl(\sum_i a_i\Bigr)! \frac{\max_i a_i}{\prod_i a_i!}. $$ Indeed, we may assume without loss of generality that the $N$-th track is heard last; conditional on this assumption, the probability that the $i$-th track will be heard $a_i$ times for each $i<N$ is $N^{-\!\sum_i a_i}$$N^{-\!\sum_i \! a_i}$ times the multinomial coefficient $(\sum_i a_i)! / \prod_i a_i!$$(\sum_i a_i)! \big/ \prod_i a_i!$. Numerically, these values are $$ 2.00000, \quad 2.84610+, \quad 3.49914-, \quad 4.02595\!- $$ for $N=2,3,4,5$.
@ A closed form for $M_N$ is available for $N \leq 3$ and probably not beyond. Trivially $M_1 = 1$; and N.Dupont already obtained the value $M_2 = 2$ by evaluating $M_2 = \sum_{a \geq 1} a/2^a$. But for $N=1$ and $N=2$ the problem reduces to the classical coupon collector's problem. Already for $N=3$ we have a surprise: $M_3 = 3/2 + (3/\sqrt{5})$, which has an elementary but somewhat tricky proof. For $N=4$, I get $$ M_4 = \frac73 - \sqrt{3} + \frac4\pi \int_{x_0}^\infty \frac{(2x+1) \, dx}{(x-1) \sqrt{4x^3-(4x-1)^2}} $$ where $x_0 = 3.43968\!+$ is the largest of the roots (all real) of the cubic $4x^3-(4x-1)^2$. I don't expect this to simplify further: itthe integral is the period over an elliptic curve of a differential with two simple poles that do not differ by a torsion point. In general one can reduce the $(N-1)$-fold sum to an $(N-2)$-fold one (which is one route to the value of $M_3$ and $M_4$), or to an $(N-3)$-fold integral, but probably not beyond.
try(N) = v=vector(N); while(!vecmin(v),v[random(N)+1]++); vecmax(v)
[turns out that one doesn't need to call vecmin each turn:
try(N, m,i)= v=vector(N); m=N; while(m, i=random(N)+1; v[i]++; if(v[i]==1,m--)); vecmax(v)
does the same thing in $\rho+O(1)$ operations per shuffle rather then $\rho+O(N)$, where $\rho$ is the cost of one random(N) call.]
So for example
averages 100001000 samples for $M_{100}$; this takes a few seconds, and seems to give about $11.7$.