  We obtain a full asymptotic expansion for orthogonal polynomials with respect
to weighted area measure on a Jordan domain $\mathscr{D}$ with real-analytic
boundary. The weight is fixed and assumed to be real-analytically smooth and
strictly positive, and for any given precision $\varkappa$, the expansion holds
with an $\mathrm{O}(N^{-\varkappa-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. In fact, 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 $\partial\mathscr{D}$.
