I hope I understood the OP correctly, in case not please let me know. And I will discuss the case of eigenvalue instead of singular value without much loss of generality.
In case you are only interested in the distribution of eigenvalues of an i.i.d. sample covariance matrix(which is also a random matrix but underlying distribution is always of dimension one), the correct reference should be [Anderson] instead of [Johnstone]. [Johnstone] dealt with a much more general case where the underlying distribution is of multi-dimension and especially with double-Wishart situation and its associated concentration bound. On the other hand [Anderson] did not consider $n(p),m(p)$ but only $m,n$ as separated parameters that tend to $\infty$, which is not suitable in "machine learning era".
$\blacksquare$1.Theorem 5.31 in [Vershynin] and Theorem 5.32 as corollary.
Unlike what the OP claimed, Bai's result also requires finite fourth moment [link1] which is not quite easy to verify in some applications, say the censored data. Assuming for the moment that the fourth moment of the underlying distribution is finite and consider the simple Gaussian regression for example, we usually study $N\approx n$ where $N$ is the sample size and $n$ is the number of variables in the model. You are right on the point if $N\gg n$ then all eigenvalues concentrate on a certain neighborhood of $\sqrt{N}$. However, a challenge in the field is that when $n$ is large and $N\approx n$. Intuitively, when the dimension of the regression model is high and the data is sparse. In this situation, some "machine learning" methods were proposed to select variables [Hastie&Tibshirani]. Among these methods, one is called sparse PCA which penalized sparse input and "forces the eigenvalues to deviate from $N$ with positive probability" if you want to view in Bai's law (Theorem 5.31).
So it is not that
"so it seems reasonable to infer that the (random) matrices people encounter in
machine learning don't behave like they have i.i.d. entries."
instead researchers will treat i.i.d. in companion with consideration of the dimension of the space that the working model lied in.
Dependent data will look almost like i.i.d data if the dimension is high. [Pastur et.al]
The consideration of dimension is not unusual even in classical statistic, some admissible models will become inadmissible if the dimension increase, for example Stein's effect [link2].
$\blacksquare$2. What kind of distribution on random matrices has this kind of distribution of singular values(eigenvalues)?
2.1 Behavior of eigenvalues of dependent Gaussian random matrices
The independent assumption is not necessary as you already noted. However, even if we dropped the independence we can still obtain a concentration bound for the largest eigenvalue. Consider $N(\boldsymbol{0},\Sigma)$ instead of $N(\boldsymbol{0},I)$. Again we consider regression setting, the following result tells you that without independence between variables we can still impose a bound on the covariance matrix which is required to proceed PCA or other variable selection techniques using eigenvalue decomposition. Notice that $Wishart(I,\bullet)$ is not a serious constraint as commented by [Johnstone] p.3.
Thm 1 in [Johnstone](modified)
Assume $p$ is even and that $p, m = m(p)$ and $n= n(p)\rightarrow\infty$ together in such a way that $\lim_{p\rightarrow\infty}\frac{\min(n,p)}{(m+n)}>0,\lim_{p\rightarrow\infty}\frac{p}{m}<1$.
For each $s_{0}\in\mathbb{R}$, there exists $C > 0$ such that for $s\geq s_0$, $$|Pr\{W_p\leq\mu_p+\sigma_p\cdot s\}-F_1(s)|<Cp^{-2/3}e^{-s/2}$$ where
$W_p=\log\left(\frac{\theta_{1,p}}{1-\theta_{1,p}}\right)$ and $\theta_{1,p}$ is
the largest eigenvalue of $(A + B)^{-1}B$ where
$A\sim Wishart_p(I,m),B\sim Wishart_p(I,n)$ are two random matrices. Here $C,\mu_p,\sigma_p$
depends ONLY on $(\gamma, \phi)$ defined as canonical angle parameters
$$\sin^2(\frac{\gamma}{2})=\frac{\min(p,n)-.5}{m+n-1}$$
$$\sin^2(\frac{\phi}{2})=\frac{\max(p,n)-.5}{m+n-1}$$ and also on
$s_0$ if $s_0 < 0$.
$F_1$ is the Tracy-Widom distribution.
To see more clearly how this will work in a regression case, we can actually see from definition that $\theta_{1,p}$ defined above is the solution to the generalized eigenvalue problem for random matrices $A,B$ such that $\det[x(A + B) − A] = 0$. A special case is $\det[P_{X}x − I] = 0$ where $P_X$ is chosen to be a projection matrix(OR the estimated covariance matrix) in regression setting. But as expected as we saw in the result above, as dimension of the model $p$ increase, such a convergence will slow down at the rate of $O(p^{-2/3})$, which allows more freedom for "extremely large/small eigen values".
I plotted an illustration as below (NOT Exact simulation!)
So for most Gaussian data, the largest eigenvalues will exponentially converge to a certain distribution, which means that the decaying rate of the largest eigenvalues is determined by the family of Tracy-Widom distribution indexed by (sample size) $m(p),n(p)$.
Moreover we can even study the speed of convergence in a more precise way using
the "central limit theorem"(circular law) for random matrices such as [Jiang].
[Bandeira&van Handel] now gives a sharp bound on spectral norm.
2.2 Behavior of eigenvalues of product uniform random matrices
Results of similar spirit concerning non-Gaussian data could be found in [Absil et.al]. In that paper the authors concern about the uniform distribution over a Grassmanian and their Theorem 1 (not copied here) also showed that the distribution of canonical angles between the spaces spanned by columns of two uniformly distributed random matrices $A,B$ are also dominated by Tracy-Widom distribution (multiplied by a gamma-polynomial term in $p,n(p)$). The cosines of canonical angles are exactly the singular values of $A'B$.[link3, Absil et.al]. An earlier related result is obtained in [Diaconis&Shahshahani] which considers only unitary random matrices.
Therefore the relative magnitude $\frac{\lambda_1}{\lambda_2}$ instead of $|\lambda_1-\lambda_2|$ are used as criterion when we choose which $k$ principle components are to remain in the model. This also explains a well-known drawback of PCA, which is that PCA is scale-dependent.
2.3 Conclusion
The asymptotic distribution of eigenvalues of a (OR asymptotically w.r.t. $m,n$) Gaussian random matrices is dominated by the Tracy-Widom distribution and its convergence behavior is controlled by an exponential rate associated with deviation $s$ and a polynomial rate associated with dimension $p$. Therefore whether the eigenvalues of a given random matrix concentrates together depends on the $m(p),n(p),p$ and eventually on the dimension $p$.
$\bullet$ The eigenvalues of a given Gaussian random matrix will not decay rapidly if we make the growth of the sample size satisfy $m(p),n(p)\rightarrow\infty$ in such a way that $$\lim_{p\rightarrow\infty}\frac{\min(n,p)}{(m+n)}>0,\lim_{p\rightarrow\infty}\frac{p}{m}<1$$. If it decays rapidly, we have some existing methods to penalized it as mentioned in part 1.
$\bullet$ The eigenvalues of product uniform distributions (on Grassmanians of dim $p$) will not decay rapidly under certain constraints.
Therefore it is safe to say that in most applications, as long as $m(p),n(p)$(usually sample sizes) grows with $p$(usually model
dimension) in a reasonable manner, the eigenvalues will neither
concentrate severely NOR change their relative magnitude too
dramatically.
On the other hand, if a random matrix is actually generated from a
model with varying (increasing) sample size $N$, then it will most
likely have a set of rapidly decaying eigenvalues under mild
constraints as discussed above.
In case where the growth conditions are not satisfied, the word "sparse" come into play.
$\blacksquare$3. Return to OP's last question.
Therefore the focus of (statistics) research is not on which family of random models have rapidly decaying eigenvalues such that we can avoid using variable selection on these "morbid models"; the focus of (some) statistical methods are on how to wrap up these possible morbidity and make variable selection based on such eigenvalues still work.
Later I have a nice discussion with some expert in the field, who pointed out that the eigenvalue problem can be framed as the distribution of roots of a polynomial with random coefficients(coefficients as random variables) and simultaneously a kernel if the covariance matrix is p.d.
The distribution of eigenvalues of a random matrix is described in terms of the distribution of roots of a determinant equation, which turns out to be a random polynomial[Deift].
And by now I think your question that
Have random models with this kind of behavior been studied, and do
their singular values behave as expected?
is addressed.
Reference
[Anderson]Anderson, Theodore Wilbur. "Asymptotic theory for principal component analysis." The Annals of Mathematical Statistics 34.1 (1963): 122-148.
[link1]https://math.stackexchange.com/questions/1467581/what-does-it-mean-for-a-probability-distribution-to-have-a-finite-fourth-moment
[Vershynin]http://www-personal.umich.edu/~Vershyninv/papers/non-asymptotic-rmt-plain.pdf
[Hastie&Tibshirani]Trevor Hastie and Robert Tibshirani "Statistical Learning with Sparsity: The Lasso and Generalizations."
[link2]https://stats.stackexchange.com/questions/263879/is-an-admissible-minimax-rule-unique
[Johnstone]Johnstone, Iain M. "Multivariate analysis and Jacobi ensembles: Largest eigenvalue, Tracy–Widom limits and rates of convergence." Annals of statistics 36.6 (2008): 2638.(Slides: http://www.stat.harvard.edu/NESS10/Johnstone-slides.pdf)
[Tracy&Widom]Tracy, Craig A., and Harold Widom. "On orthogonal and symplectic matrix ensembles." Communications in Mathematical Physics 177.3 (1996): 727-754.
[Absil et.al]Absil, P-A., Alan Edelman, and Plamen Koev. "On the largest principal angle between random subspaces." Linear Algebra and its applications 414.1 (2006): 288-294.
[link3]angle between subspaces
[Jiang]Jiang, Tiefeng. "Approximation of Haar distributed matrices and limiting distributions of eigenvalues of Jacobi ensembles." Probability theory and related fields 144.1-2 (2009): 221-246.
[Diaconis&Shahshahani]Diaconis, Persi, and Mehrdad Shahshahani. "On the eigenvalues of random matrices." Journal of Applied Probability (1994): 49-62.http://statweb.stanford.edu/~cgates/PERSI/papers/random_matrices.pdf
[Pastur et.al]Pastur, Leonid Andreevich, Mariya Shcherbina, and Mariya Shcherbina. Eigenvalue distribution of large random matrices. Vol. 171. Providence, RI: American Mathematical Society, 2011.
[Deift]Deift, Percy. Orthogonal polynomials and random matrices: a
Riemann-Hilbert approach. Vol. 3. American Mathematical Soc., 2000.
[Bandeira&van Handel]Bandeira, Afonso S., and Ramon van Handel. "Sharp nonasymptotic bounds on the norm of random matrices with independent entries." The Annals of Probability 44.4 (2016): 2479-2506.
Best Answer
The "Katz-Sarnak philosophy" is just the idea that statistics of various kinds for $L$-functions should, in the large scale limit, match statistics for large random matrices from some particular classical compact group.
First you need to decide what kinds of zeros to look at: the high zeros of an individual $L$-function, the low zeros (near the real axis) in some family of related $L$-functions (e.g., low zeros in all Dirichlet $L$-functions, maybe with a restriction on the parity of the character), and so on. The choice of what sort of zeros you look at can have an effect, e.g., looking at statistics of the low zeros in a family can reveal a "symmetry" (similar statistics to a specific classical compact group) that is not evident when looking at high zeros of an individual $L$-function (universality, not distinguishing between different $L$-functions).
Next you need to normalize the (nontrivial) zeros. For example, zeros of an individual $L$-function tend to get close to each other high up the critical line, so you count the number of zeros up to height $T$ and then rescale them so they get average spacing 1. It's those rescaled zeros that you work with when computing your various statistics.
Next you need to work with a suitable class of test functions and be able to actually carry out statistical calculations to reveal similarity with some class of random matrices.
Part of the motivation for people to look at these questions is the hope that it could suggest a spectral interpretation of the nontrivial zeros. Look up the Hilbert-Polya conjecture (Hilbert had nothing to do with it).
You mentioned the Katz-Sarnak book in your question. Keep in mind that this book is entirely about $L$-functions for varieties over finite fields, not the ones like $\zeta(s)$ that are associated to number fields. The blog post and paper you mention by Kowalski and Miller are about $L$-functions over number fields, and that is perhaps why they look very different to you from the Katz-Sarnak book. Another paper where the number field case is treated head-on is Rudnick, Sarnak, "Zeros of principal $L$-functions and random matrix theory," Duke Math. J. 81 (1996), 269–322. If you want to see a textbook's treatment of the relation between $L$-functions over number fields (like $\zeta(s)$ and Dirichlet $L$-functions) and random matrix theory, consider the book "An Invitation to Modern Number Theory" by Steve Miller and Ramin Takloo-Bighash. The random matrix theory is in Part 5.