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?

See also this application.

How I got this far

first term

For the first term, start with the integral representation

$$ 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$:



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) -

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.

  • 1
    $\begingroup$ Nitpick: this is not a Laurent series, it's a Puiseux series $\endgroup$
    – Wojowu
    Mar 26 '19 at 8:46
  • $\begingroup$ Thanks @Wojowu, edited & learned something. $\endgroup$
    – Roman
    Mar 26 '19 at 8:49
  • $\begingroup$ I don't know what has been your process (I really would like to know) but this is beautiful. $\endgroup$ Mar 26 '19 at 9:07
  • 1
    $\begingroup$ Have you seen this? $\endgroup$ Mar 26 '19 at 11:25
  • 1
    $\begingroup$ Yes, Watson indeed has it; the formula is attributed to Meissel. $\endgroup$ Mar 26 '19 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

λ[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_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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.