  Up to now, the nonparametric analysis of multidimensional continuous-time
Markov processes has focussed strongly on specific model choices, mostly
related to symmetry of the semigroup. While this approach allows to study the
performance of estimators for the characteristics of the process in the minimax
sense, it restricts the applicability of results to a rather constrained set of
stochastic processes and in particular hardly allows incorporating jump
structures. As a consequence, for many models of applied and theoretical
interest, no statement can be made about the robustness of typical statistical
procedures beyond the beautiful, but limited framework available in the
literature. To close this gap, we identify $\beta$-mixing of the process and
heat kernel bounds on the transition density as a suitable combination to
obtain $\sup$-norm and $L^2$ kernel invariant density estimation rates matching
the case of reversible multidimenisonal diffusion processes and outperforming
density estimation based on discrete i.i.d. or weakly dependent data. Moreover,
we demonstrate how up to $\log$-terms, optimal $\sup$-norm adaptive invariant
density estimation can be achieved within our general framework based on tight
uniform moment bounds and deviation inequalities for empirical processes
associated to additive functionals of Markov processes. The underlying
assumptions are verifiable with classical tools from stability theory of
continuous time Markov processes and PDE techniques, which opens the door to
evaluate statistical performance for a vast amount of Markov models. We
highlight this point by showing how multidimensional jump SDEs with L\'evy
driven jump part under different coefficient assumptions can be seamlessly
integrated into our framework, thus establishing novel adaptive $\sup$-norm
estimation rates for this class of processes.
