As per the title: consider an $n\times n$ matrix of the form $$J_n = a I_n + E_n,$$ where $E_n$ is defined componentwise as $(E_n)_{i,i+1}=1$ for all $i$, and with zero components everywhere else. This is a single Jordan block of size $n$.
Is there a general expression for its singular values?
This is simple enough to work out in the $n=2$ case, where we get $$s_\pm(J_2) = \frac{1}{\sqrt2}\sqrt{1 + 2|a|^2 \pm \sqrt{1 + 4 |a|^2}}.$$ Already for $n=3$ I get nontrivial cubic polynomial equations to solve, however. While of course it is natural that working this out for generic $n$ involves solving polynomial equations of degree $n$, is there any way to simplify the problem, or get some understanding about the structure of the solutions more in general?
Looking at it numerically, there seems to be a relatively straightforward structure to the solutions. For example, for $n=4$ and $a\in\mathbb{R}$, we have
where I'm plotting the singular values of $J_4$ as a function of $a\in[-3,3]$. It also seems like the only thing that ever matters is $|a|$, so the right half-plane in this plot is sufficient to cover the general case with $n=4$, also when $a\in\mathbb{C}$. There is clearly a lot of structure here, so how can we see it analytically by computing or somehow estimating the singular values?

I'll write the nilpotent part of the Jordan block as $N$. You can compute that $\| N \| = 1$, so $J$ is within $1$ in the operator norm of the scalar $aI$, from which it follows (because the singular values are Lipschitz with Lipschitz constant $1$) that the singular values of $J$ are within $1$ of $|a|$, which explains the rough shape of the asymptotic behavior as $|a| \to \infty$.
To get something more specific, $J^{\dagger} J$ is almost the matrix $M = (|a|^2 + 1) I + a N^{\dagger} + \bar{a} N$, except that the top left entry is $|a|^2$ instead of $|a|^2 + 1$. $M$ is a tridiagonal Toeplitz matrix, which has eigenvalues
$$\lambda_i = |a|^2 + 2 |a| \cos \frac{i \pi}{n+1} + 1, 1 \le i \le n$$
(the eigenvectors can also be written down explicitly in terms of trigonometric functions). $J^{\dagger} J$ is a perturbation of this matrix by a matrix of operator norm $1$ as above, so we get (using the Lipschitz property of singular values again, and using that $J^{\dagger} J$ is positive-semidefinite and Hermitian, meaning its eigenvalues and singular values agree) $|\sigma_i^2 - \lambda_i| \le 1$, which gives that the singular values of $J$ have asymptotic behavior
$$\sigma_i = |a| + \cos \frac{i \pi}{n+1} + O \left( \frac{1}{|a|} \right)$$
as $|a| \to \infty$. To match this up with your computations, when $n = 4$ this gives singular values which are roughly $|a|$ plus, respectively,
$$\cos \frac{\pi}{5} \approx 0.809 \dots $$ $$\cos \frac{2\pi}{5} \approx 0.309 \dots $$ $$\cos \frac{3 \pi}{5} \approx -0.309 \dots $$ $$\cos \frac{4 \pi}{5} \approx -0.809 \dots $$
which roughly agrees visually with your plot at the edges.