  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}$.
