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 β ensemble, characterized by the Dyson parameter β, and the Laguerre weight xae−β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↦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∈double-struckZ≥0. We show, using the theory of multidimensional hypergeometric functions based on Jack polynomials, that with the modified hard edge scaling x↦x/4(N+a/β), the rate of convergence to the limiting distribution is O(1/N2), which is optimal. In the case β=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 β>0. An iterative scheme is presented to numerically approximate the functional form for general a∈double-struckZ≥2.