What is the Puiseux series of the Bessel function $J_n(n)$?

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

$$J_n(n) = \frac{\Gamma(\frac13)}{2^{2/3}\cdot 3^{1/6} \cdot \pi}n^{-1/3} - \frac{1}{35\cdot 6^{1/3}\cdot\Gamma(\frac13)}n^{-5/3} - \frac{\Gamma(\frac13)}{225 \cdot 2^{2/3}\cdot 3^{1/6}\cdot \pi}n^{-7/3} +\mathcal{O}(n^{-11/3})$$

How does this series continue?

How I got this far

first term

$$J_n(n) = \frac{1}{2\pi}\int_{-\pi}^{\pi} d\theta \cos[n(\sin(\theta)-\theta)]$$

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

$$\lim_{n\to\infty} n^{1/3}\cdot\frac{1}{2\pi}\int_{-\pi}^{\pi} d\theta \cos[-n\theta^3/6] = \frac{\Gamma(\frac13)}{2^{2/3}\cdot3^{1/6}\cdot\pi}$$

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 $$n\to\infty$$:

NumericalMathNSequenceLimit[T]


-0.00586928848357833870

Then use AskConstants to find that this number is probably equal to $$-\frac{1}{35\cdot 6^{1/3}\cdot\Gamma(\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 $$n\to\infty$$. The result is $$-0.0019880325262065435671$$, which AskConstants thinks is equal to $$- \frac{\Gamma(\frac13)}{225 \cdot 2^{2/3}\cdot 3^{1/6}\cdot \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

$$+\frac{1213}{511875\cdot6^{1/3}\cdot \Gamma(\frac13)}n^{-11/3}$$

but such large rationals don't inspire confidence.

• Nitpick: this is not a Laurent series, it's a Puiseux series – Wojowu Mar 26 at 8:46
• Thanks @Wojowu, edited & learned something. – Roman Mar 26 at 8:49
• I don't know what has been your process (I really would like to know) but this is beautiful. – Claude Leibovici Mar 26 at 9:07
• Have you seen this? – J. M. is a poor mathematician Mar 26 at 11:25
• Yes, Watson indeed has it; the formula is attributed to Meissel. – J. M. is a poor mathematician Mar 26 at 11:51

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

$$J_n(n) = \frac{1}{\pi} \sum_{m=0}^{\infty} \lambda_m \cdot \Gamma\left(\frac{2m+4}{3}\right) \cdot \left(\frac{6}{n}\right)^{\frac{2m+1}{3}} \cdot \cos\left(\frac{2m+1}{6}\pi\right)$$

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=\frac{1}{60}\\ \lambda_2=\frac{1}{1400}\\ \lambda_3=\frac{1}{25200}\\ \lambda_4=\frac{43}{17248000}\\ \lambda_5=\frac{1213}{7207200000}\\ \lambda_6=\frac{151439}{12713500800000}\\ \lambda_7=\frac{33227}{38118080000000}\\ \lambda_8=\frac{16542537833}{252957982717440000000}\\ \lambda_9=\frac{887278009}{177399104762880000000}\\ \lambda_{10}=\frac{15233801224559}{39217856135377920000000000}\\ \ldots$$

These coefficients can be calculated efficiently with the Mathematica code

λ = 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_0\ldots\lambda_n$$ with

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


$$x+\frac{x^3}{60}+\frac{x^5}{1400}+\frac{x^7}{25200}+\frac{43 x^9}{17248000}+\frac{1213x^{11}}{7207200000}+\mathcal{O}(x^{12})$$

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

• If you need a whole pile at once: ComposeSeries[InverseSeries[Series[u - Sin[u], {u, 0, 25}]], x^3/6 + O[x]^25] – J. M. is a poor mathematician Mar 26 at 14:19
• Thanks @J.M.isnotamathematician , I was wondering how to make InverseSeries` work on this case. – Roman Mar 26 at 14:35