Let $c(n, k)$ denote the number of permutations in $S_n$ whose cycle decomposition has $k$ cycles. For a fixed $n$, I want to find $k$ such that $c(n, k)$ is maximized.
I know that the $k$ I seek is either the floor or ceiling of $1+\frac{1}{2}+\cdots+\frac{1}{n}$, but I'm not quite sure how to do this. This maybe suggests using an exponential generating function, but I don't know what the EGF of the sequence $c(n, k)$, for a fixed $n$ is. Could I have some hints?
I’m not sure how useful this is for you, but in this paper [PDF] Erdős showed that there is a $k_n$ maximizing $c(n,k)$ and that it is
$$k_n=\left\lfloor\ln(n+1)+\gamma-1+\frac{\zeta(2)-\zeta(3)}{\ln(n+1)+\gamma-\frac32}+\frac{h}{\left(\ln(n+1)+\gamma-\frac32\right)^2}\right\rfloor\;,$$
where $\gamma$ is the Euler-Mascheroni constant, $\zeta$ is the Riemann zeta function, and $-1.1<h<1.5$. He also noted that for $n>188$ this can be simplified to
$$\left\lfloor\ln n-\frac12\right\rfloor\le k_n\le\lfloor\ln n\rfloor\;.$$
I have not looked at the Erdös paper yet but I would like to thank Brian M. Scott for the quick reply giving such a precise reference.
What follows are some informal observations where the reader is asked to be patient with the lack of rigor. Some time ago I used the Polya Enumeration Theorem to prove the following asymptotic for Stirling numbers of the first kind (this is the MSE link) $$\left[ n\atop k\right] \sim \frac{(n-1)!}{(k-1)!} \log^{k-1} (n-1).$$
We ask for what $k$ this is maximized. Consider the function $$\frac{Q^x}{\Gamma(x+1)}$$
with $Q>1$ a constant . We have by inspection that the growth from the exponential term dominates until the Gamma function term takes over, for an ultimate limit of zero.
To locate the point where this happens differentiate to get $$\log Q \times \frac{Q^x}{\Gamma(x+1)} - \frac{Q^x}{\Gamma(x+1)^2} \Gamma'(x+1) = 0.$$ This gives $$\log Q - \frac{\Gamma'(x+1)}{\Gamma(x+1)} = 0$$ or in terms of the digamma function $$\psi(x+1) = \log Q.$$ But on the real line we have $\psi(x) \sim \log x$ so that the conclusion is that $$x\sim Q.$$
Returning to the Stirling numbers we see that here $Q = \log(n-1)$ giving the approximation $$k\sim\log n$$ for the $k$ that maximizes $\left[n\atop k\right]$ with $n$ fixed.
Addendum. Where generating functions are concerned we have the species $$\mathfrak{P}(\mathcal{U}\mathfrak{C}(\mathcal{Z}))$$ which gives the bivariate generating function $$\exp\left(u\log\frac{1}{1-z}\right)$$ so that with $n$ fixed $$\left[n\atop k\right] = n! [z^n] \frac{1}{k!} \left(\log\frac{1}{1-z}\right)^k.$$ The mechanics of extracting coefficient asymptotics from this are discussed in the text Analytic Combinatorics by Flajolet and Sedgewick and in the slides from that text which refer to the so-called standard function scale.
Remark. The fomulae from the PET computation are exact, even though we have used only the first term here.
