Using
$$\int J_0(x)\,dx =x \, _1F_2\left(\frac{1}{2};1,\frac{3}{2};-\frac{x^2}{4}\right)$$
$$\int_0^a J_0(x)\,dx =\frac{1}{2} a (\pi \pmb{H}_0(a) J_1(a)+(2-\pi \pmb{H}_1(a)) J_0(a))$$
$$\int_0^{j_{0,n}} J_0(x)\,dx=\frac{\pi}{2}\, j_{0,n}\, \pmb{H}_0\left(j_{0,n}\right)\, J_1\left(j_{0,n}\right)$$
$$\sigma_n=(1)^n\frac{\pi}{2} \big( j_{0,n+1} \,\pmb{H}_0\left(j_{0,n+1}\right)\,
J_1\left(j_{0,n+1}\right)- j_{0,n}\, \pmb{H}_0\left(j_{0,n}\right)\,
J_1\left(j_{0,n}\right)\big)$$
it seems to be that
$$\color{blue}{\sigma_n\,=\frac{8}{\pi\sqrt {8n}}\left(1-\frac{1}{(8 n)}+\frac{1}{(8n)^2} -\frac{1}{10}\frac{1}{(8n)^3}-\frac{2}{(8n)^4}+O\!\left(\frac1{n^5}\right)\right)}$$
This was obtained for values up to $n=100$.
For a few small values of $n$, here are the results
$$\left(
\begin{array}{ccc}
n & \text{approximation} & \text{exact} \\
1 & 0.8012287685 & 0.8014542111 \\
2 & 0.5992828620 & 0.5993225154 \\
3 & 0.4990351563 & 0.4990496204 \\
4 & 0.4365280908 & 0.4365351123 \\
5 & 0.3928185568 & 0.3928225593 \\
6 & 0.3600543084 & 0.3600568365 \\
7 & 0.3343192651 & 0.3343209797 \\
8 & 0.3134138472 & 0.3134150722 \\
9 & 0.2959950956 & 0.2959960063 \\
10 & 0.2811906203 & 0.2811913190 \\
20 & 0.2000664765 & 0.2000665991 \\
30 & 0.1636924770 & 0.1636925214 \\
40 & 0.1419090468 & 0.1419090684 \\
50 & 0.1270064402 & 0.1270064525 \\
60 & 0.1159886945 & 0.1159887023 \\
70 & 0.1074165671 & 0.1074165724 \\
80 & 0.1005013911 & 0.1005013949 \\
90 & 0.0947700476 & 0.0947700505 \\
100 & 0.0899192327 & 0.0899192349
\end{array}
\right)$$
Edit
Thanks to @Roman's answer to this question of mine, we now have the coefficients for
$$\sigma_n=\frac{2\sqrt2}{\pi\sqrt n}\left(1+\sum_{k=1}^{18}\frac{a_k}{n^k}+O\left(\frac1{n^{19}}\right)\right)$$
Using these coefficients, the errors are smaller than $10^{-10}$ as soon as $n >3$. For $n=10$, the error is $4 \times 10^{-20}$ and for $n=100$, the error is $3 \times 10^{-29}$.
The following might help, but it's not a complete solution.
The coefficients $a_n$ have an alternative expression, where the $B_n$ are the Bernoulli numbers:
$$ \tag{1} a_n = \sum_{k=0}^n 2^k a_k \ \frac{B_{n-k}}{(n-k)!} $$
To prove $(1)$, start with the recursion given in the problem statement, flip the direction of summation, and pull the $a_n$ term into the sum to get
$$2^n a_n = \sum_{k=0}^n \frac{1}{k+1} \binom{n}{k} a_{n-k} \frac{(n-k)!}{n!} $$
Set $c_{n} = n!a_n$ so that
$$2^n a_n = \sum_{k=0}^n \frac{1}{k+1} \binom{n}{k} c_{n-k} $$
This equation is invertible, e.g. Binomial Transform of a Particular Binomial Sum
Reverting back to $a_n$ gives the formula $(1)$. Define a generating function $f(x) =\sum_{n \ge 0} a_n x^n$. Then $(1)$ implies the functional equation
$$ \tag{2} \quad f(x) = f(2x) \frac{x}{e^x -1} $$
Iterate on $(2)$, and formally the solution is
$$ \tag{3} \quad f(x) = \prod_{k=0}^\infty \frac{ \exp{\big(x \ 2^{-k-1}\big)} - 1}{x \ 2^{-k-1} } $$
I've calculated many of the $a_n$ and checked that $(3)$ is satisfied numerically when the product is truncated and $x$ is sufficiently small.
$(3)$ looks like a formidable equation to do asymptotic analysis on. $(2)$ looks like a trick with a Mellin transform might be advantageous. The closest thing I know of that looks like the proposed asymptotic expansion has to do with Rice integrals, and instead of Bernoulli numbers, binomials appear in the sum. See example $4$ in $[1]$.
Edit 1/13/2022.
Summary: the asymptotics of $a_n$ can be approximated by a saddle point analysis of the large $x$ asymptotics of $f(x)$ given above. Here is a short algorithm for calculation:
$$ \tag{4} x_0=-\frac{1}{\log{2}}W_{-1}\big( -\log{2} \cdot 2^{-n} \big) $$
$$ \tag{5} C=2^{-1/12}\exp{\Big(\big( \gamma_1 + \gamma^2/2-\pi^2/12 \big) }/\log{2} \Big) $$
$$ \tag{6} F(x)=C\sqrt{x}\exp{\big(x -\frac{1}{2 \log{2}}\log^2 x \big)}$$
$$b(x)=x \frac{d}{dx}(x \frac{d}{dx} \log{(F(x)} ) $$
$$\tag{7} a_n^{A} \sim (2 \pi b(x_0))^{-1/2} F(x_0) x_0^{-n} $$
Notation: $W_{-1}(x)$ is the Lambert W function on the '-1 branch' and is implemented in Mathematica as ProductLog[-1,x]. The $\gamma$ is Euler's constant and $\gamma_1$ is the first Stieltjes constant, $\gamma_1 \approx$ -0.0728158.
Here is a table comparing the true $a_n,$ the approximation (above) $a_n^{A},$ and the proposed approximation, $a_n^{P} = 2^{ (eq.3)}$ where eq. 3 corresponds to that in the proposers write up.
$$
\begin{array}{c|lcr}
n & a_n & a_n^{A} & a_n^{P} \quad \\
\hline
20 & 6.08905 \,\cdot\, 10^{-22} & 6.07191\,\cdot\, 10^{-22} & 9.3967\,\cdot\, 10^{-22} \\
40 & 1.00137 \,\cdot\, 10^{-52}& 1.00047\,\cdot\, 10^{-52} & 1.3980\,\cdot\, 10^{-52} \\
60 & 1.30647\,\cdot\, 10^{-87} & 1.30565\,\cdot\, 10^{-87} & 1.7222\,\cdot\, 10^{-87} \\
80 & 3.13549\,\cdot\, 10^{-125} & 3.13398\,\cdot\, 10^{-125} & 3.9804\,\cdot\, 10^{-125}
\end{array}
$$
The approximation $a_n^{A}$ is excellent, about 3 significant digits for $n=20,$ and about 4 digits for $n=80.$ A comment before sketching the proof. The proposer's solution seems to have missed the constant, and otherwise the expansion seems to be an asymptotic expansion of the Lambert W function,
$$ W(x) \sim \log{x} - \log{\log{x}} + \cdots $$
The first part of the proof is to show that, for large $x$
$$ \tag{4} \quad f(x) = \prod_{k=1}^\infty \frac{ \exp{\big(x \ 2^{-k}\big)} - 1}{x \ 2^{-k} }
\sim F(x)\exp{\big(P(x) \big)}
$$
where $F(x)$ is given by eqs. (5) & (6) and $P(x)$ is given by
$$P(x)=\sum_{k=-\infty, k \ne 0}^\infty\ \frac{1}{2\pi i k}\Gamma\big(1-\frac{2 \pi i k}{\log2} \big)
\zeta\big(1-\frac{2 \pi i k}{\log2} \big) \exp{\big(x \, \frac{2 \pi i k}{\log2} \big)}.$$
The formulas presented are a slight modification of that presented in Example 9.3 in reference [2]. We are going to ignore $P(x),$ in subsequent analysis, mainly because the amplitudes of the oscillations are tiny, with magnitudes of $7.4\,\cdot\,10^{-7}, 4.4\,\cdot\,10^{-13}\cdots$ for $k=1,2...$ However, if you want even more precise asymptotics (3 or 4 digits is not enough?) then they have to be included.
Reference [2] deals with a problem that has the functional equation
$$ g(x) = g(2x)\frac{x}{1-e^{-x}} $$
Note how close that is to the functional equation in (3). On taking logarithms we have
$$\log{\big(f(x)\big)} = \sum_{k=1}^\infty \log{\Big(
\frac{ \exp{\big(x \ 2^{-k}\big)} - 1}{x \ 2^{-k} } \Big) }
$$
$$ = \sum_{k=1}^\infty \log{\Big(\exp{\big(x \ 2^{-k}\big)}\Big) } + \log{\Big(
\frac{ 1-\exp{\big(-x \ 2^{-k}\big)} }{x \ 2^{-k} } \Big) }
$$
$$=x+ \sum_{k=1}^\infty \log{\Big(
\frac{ 1-\exp{\big(-x \ 2^{-k}\big)} }{x \ 2^{-k} } \Big) }$$
The infinite sum is precisely that given in ref. [2] for $g(x)$, so its asymptotics can be copied verbatim, which I have done. Very briefly, the technique is to find the Mellin transform for the sum, then find the asymptotics by moving the contour of integration over the function's singularities. Later on, the reference extracts the asymptotics for the coefficients of $g(x)$ by using a process called dePoissionization, but that requires a growth restriction on the function. The addition of $e^x$ for my $f(x)$ destroys that possibility. However, traditional saddle point works. In that situation, one can use eq. (7) for a generic function $F(x).$ How does one determine the saddle point? Look at the Cauchy representation of $a_n,$
$$ a_n = \frac{1}{2 \pi i} \oint \frac{F(x)}{x^{n+1}} dx$$
$$ \rightarrow \frac{1}{2 \pi i} \oint \exp{\big( x - n\log{x} - \frac{1}{2 \log{2}}\log^2{x} \big)} \frac{dx}{\sqrt{x}}
$$
The saddle point $x0$ is that in which the derivative of the argument of the exponential is zero,
$$ \frac{d}{dx} \big(x - n\log{x} - \frac{1}{2 \log{2}}\log^2{x} \big)\overset{!}{=} 0 $$
The solution to this equation for $x=x_0$ is given in eq. (4).
$[1]$ Mellin transform and asymptotics: Finite differences and Rice's integrals, P. Flajolet & R. Sedgewick, Ther. Comp. Sci 144 (1995) 101-124.
$[2]$ Average Case Analysis of Algorithms on Sequences W. Szpankowski, Wiley-Interscience, 2001.
Best Answer
Let us write $$ \sum\limits_{k = 1}^n {k^n } = n^n \sum\limits_{k = 0}^{n - 1} {\left( {1 - \frac{k}{n}} \right)^n } . $$ Now by power series expansions \begin{align*} \left( {1 - \frac{k}{n}} \right)^n & = \exp \left( {n\log \left( {1 - \frac{k}{n}} \right)} \right) = \exp \left( { - n\sum\limits_{j = 1}^\infty {\frac{{k^j }}{j}\frac{1}{{n^j }}} } \right) \\ & = \mathrm{e}^{ - k} - \mathrm{e}^{ - k} \frac{{k^2 }}{2}\frac{1}{n} + \mathrm{e}^{ - k} \left( {\frac{{k^4 }}{8} - \frac{{k^3 }}{3}} \right)\frac{1}{{n^2 }} - \ldots \,. \end{align*} Substituting back to the original sum, extending that sum to $n=\infty$, and using closed forms of $\sum_{k=0}^\infty k^px^k$ (with $x=\mathrm{e}^{-1}$), we deduce \begin{align*} \sum\limits_{k = 1}^n {k^n } & \sim n^n \left( {\frac{\mathrm{e}}{{\mathrm{e} - 1}} - \frac{{\mathrm{e}^2 + \mathrm{e}}}{{2(\mathrm{e} - 1)^3 }}\frac{1}{n} - \frac{{\mathrm{e}(5\mathrm{e}^3 - 9\mathrm{e}^2 - 57\mathrm{e} - 11)}}{{24(\mathrm{e} - 1)^5 }}\frac{1}{{n^2 }} + \ldots } \right) \\ & =\frac{\mathrm{e}}{{\mathrm{e} - 1}}n^n \left( {1 - \frac{{\mathrm{e} + 1}}{{2(\mathrm{e} - 1)^2 }}\frac{1}{n} - \frac{{5\mathrm{e}^3 - 9\mathrm{e}^2 - 57\mathrm{e} - 11}}{{24(\mathrm{e} - 1)^4 }}\frac{1}{{n^2 }} + \ldots } \right) \end{align*} as $n\to +\infty$. You can obtain more terms of the expansion if you like, the process should be clear.