We obtain a full asymptotic expansion for orthogonal polynomials with respect to weighted area measure on a Jordan domain D with real-analytic boundary. The weight is fixed and assumed to be real-analytically smooth and strictly positive, and for any given precision κ, the expansion holds with an O(N−κ−1) error in N-dependent neighborhoods of the exterior region as the degree N tends to infinity. The main ingredient is the derivation and analysis of Riemann-Hilbert hierarchies—sequences of scalar Riemann-Hilbert problems—which allows us to express all higher order correction terms in closed form. Indeed, the expansion may be understood as a Neumann series involving an explicit operator. The expansion theorem leads to a semiclassical asymptotic expansion of the corresponding hard edge probability wave function in terms of distributions supported on ∂D.