Finite size corrections at the hard edge for the Laguerre $\beta$ ensemble

A fundamental question in random matrix theory is to quantify the optimal rate of convergence to universal laws. We take up this problem for the Laguerre $\beta$ ensemble, characterised by the Dyson parameter $\beta$, and the Laguerre weight $x^a e^{-\beta x/2}$, $x > 0$ in the hard edge limit. The latter relates to the eigenvalues in the vicinity of the origin in the scaled variable $x \mapsto x/4N$. Previous work has established the corresponding functional form of various statistical quantities --- for example the distribution of the smallest eigenvalue, provided that $a \in \mathbb Z_{\ge 0}$. We show, using the theory of multidimensional hypergeometric functions based on Jack polynomials, that with the modified hard edge scaling $x \mapsto x/4(N+a/\beta)$, the rate of convergence to the limiting distribution is $O(1/N^2)$, which is optimal. In the case $\beta = 2$, general $a> -1$ the explicit functional form of the distribution of the smallest eigenvalue at this order can be computed, as it can for $a=1$ and general $\beta > 0$. An iterative scheme is presented to numerically approximate the functional form for general $a \in \mathbb Z_{\ge 2}$.