We develop a convergence analysis of a multi-level algorithm combining higher
order quasi-Monte Carlo (QMC) quadratures with general Petrov-Galerkin
discretizations of countably affine parametric operator equations of elliptic
and parabolic type, extending both the multi-level first order analysis in
[\emph{F.Y.~Kuo, Ch.~Schwab, and I.H.~Sloan, Multi-level quasi-Monte Carlo
finite element methods for a class of elliptic partial differential equations
with random coefficient} (in review)] and the single level higher order
analysis in [\emph{J.~Dick, F.Y.~Kuo, Q.T.~Le~Gia, D.~Nuyens, and Ch.~Schwab,
Higher order QMC Galerkin discretization for parametric operator equations} (in
review)]. We cover, in particular, both definite as well as indefinite,
strongly elliptic systems of partial differential equations (PDEs) in
non-smooth domains, and discuss in detail the impact of higher order
derivatives of {\KL} eigenfunctions in the parametrization of random PDE inputs
on the convergence results. Based on our \emph{a-priori} error bounds, concrete
choices of algorithm parameters are proposed in order to achieve a prescribed
accuracy under minimal computational work. Problem classes and sufficient
conditions on data are identified where multi-level higher order QMC
Petrov-Galerkin algorithms outperform the corresponding single level versions
of these algorithms. Numerical experiments confirm the theoretical results