Fdtxjr

By numerical experimentation I find the first three terms of the Puiseux series of the Bessel function of the first kind

$$
J_n(n) =
fracGamma(frac13)2^2/3cdot 3^1/6 cdot pin^-1/3
- frac135cdot 6^1/3cdotGamma(frac13)n^-5/3
- fracGamma(frac13)225 cdot 2^2/3cdot 3^1/6cdot pin^-7/3
+mathcalO(n^-11/3)
$$

How does this series continue?

See also this application.

How I got this far

first term

For the first term, start with the integral representation

$$
J_n(n) = frac12piint_-pi^pi dtheta cos[n(sin(theta)-theta)]
$$

For $ntoinfty$ the only significant contributions to this integral come from values of $theta$ that are close to zero. Therefore we approximate $sin(theta)-thetaapprox-theta^3/6$ and find

$$
lim_ntoinfty n^1/3cdotfrac12piint_-pi^pi dtheta cos[-ntheta^3/6] = fracGamma(frac13)2^2/3cdot3^1/6cdotpi
$$

In Mathematica:

Limit[1/(2π) Integrate[Cos[n (-(θ^3/6))], θ, -π, π]*n^(1/3), n -> ∞]

Gamma[1/3]/(2^(2/3) 3^(1/6) π)

second term

In Mathematica, define the Bessel function and its one-term approximation, as well as their numerical difference evaluated to 1000 digits:

b[n_] = BesselJ[n, n];
ba[n_] = Gamma[1/3]/(2^(2/3)*3^(1/6)*π)*n^(-1/3);
B[n_] := N[b[n] - ba[n], 10^3]

Calculate how the numerical difference behaves for large $n$ (after multiplying it by $n^5/3$):

ListLinePlot[T = Table[B[n]*n^(5/3), n, 10^Range[2, 5, 1/4]]]

and find the approximate numerical value of the limit as $ntoinfty$:

NumericalMath`NSequenceLimit[T]

-0.00586928848357833870

Then use AskConstants to find that this number is probably equal to $-frac135cdot 6^1/3cdotGamma(frac13)$.

third term

Same procedure as second term, but with the better approximation

ba[n_] = Gamma[1/3]/(2^(2/3)*3^(1/6)*π)*n^(-1/3) -
 1/(35*6^(1/3)*Gamma[1/3])*n^(-5/3);

and multiplying the difference B[n] by $n^7/3$ before taking the numerical limit $ntoinfty$. The result is $-0.0019880325262065435671$, which AskConstants thinks is equal to $- fracGamma(frac13)225 cdot 2^2/3cdot 3^1/6cdot pi$.

higher-order terms

The above recipe can be continued to higher-order terms, but I lose confidence in the capabilities of AskConstants. The fourth term is $+0.00048679979012516409164$, which may be

$$
+frac1213511875cdot6^1/3cdot Gamma(frac13)n^-11/3
$$

but such large rationals don't inspire confidence.

The answer is given by Ernst Meissel in this 1891 paper (in German):

$$
J_n(n) = frac1pi
sum_m=0^infty
lambda_m cdot
Gammaleft(frac2m+43right)
cdot left(frac6nright)^frac2m+13
cdot cosleft(frac2m+16piright)
$$

The coefficients $lambda_m$ describe the Taylor series of the solution $u(x)=sum_m=0^infty lambda_m x^2m+1$ of the transcendental equation $u-sin(u)=x^3/6$ around $x=0$. Term-by-term comparison gives

$$
lambda_0=1\
lambda_1=frac160\
lambda_2=frac11400\
lambda_3=frac125200\
lambda_4=frac4317248000\
lambda_5=frac12137207200000\
lambda_6=frac15143912713500800000\
lambda_7=frac3322738118080000000\
lambda_8=frac16542537833252957982717440000000\
lambda_9=frac887278009177399104762880000000\
lambda_10=frac1523380122455939217856135377920000000000\
ldots
$$

These coefficients can be calculated efficiently with the Mathematica code

λ[0] = 1;
λ[m_Integer /; m >= 1] := λ[m] = Module[Λ, u, x,
 u = Sum[λ[j] x^(2 j + 1), j, 0, m - 1] + Λ x^(2 m + 1);
 Λ /. First[Solve[SeriesCoefficient[u - Sin[u], x, 0, 2 m + 3] == 0, Λ]]]

Or all at once by series inversion (thanks to J.M.): calculate $lambda_0ldotslambda_n$ with

With[n = 5,
 ComposeSeries[InverseSeries[Series[u-Sin[u], u,0,2n+3]], x^3/6 + O[x]^(2n+5)]]

$$
x+fracx^360+fracx^51400+fracx^725200+frac43 x^917248000+frac1213x^117207200000+mathcalO(x^12)
$$

Thanks to J.M. who pointed out Meissel's paper to me.