  We study the Fredholm determinant of an integrable operator acting on the
interval $(0,s)$ whose kernel is constructed out of a hierarchy of higher order
analogues to the Painlev\'{e} III equation. This Fredholm determinant describes
the critical behavior of the eigenvalue gap probability at the hard edge of
unitary invariant random matrix ensembles perturbed by poles of order $k$ in
the double scaling regime. Using the Riemann-Hilbert method, we obtain the
large $s$ asymptotics of the Fredholm determinant. Moreover, we derive a
Painlev\'e type formula of the Fredholm determinant, which is expressed in
terms of an explicit integral involving a solution to the coupled Painlev\'e
III system.
