  We develop a large-scale regularity theory of higher order for
divergence-form elliptic equations with heterogeneous coefficient fields $a$ in
the context of stochastic homogenization. The large-scale regularity of
$a$-harmonic functions is encoded by Liouville principles: The space of
$a$-harmonic functions that grow at most like a polynomial of degree $k$ has
the same dimension as in the constant-coefficient case. This result can be seen
as the qualitative side of a large-scale $C^{k,\alpha}$-regularity theory,
which in the present work is developed in the form of a corresponding
$C^{k,\alpha}$-"excess decay" estimate: For a given $a$-harmonic function $u$
on a ball $B_R$, its energy distance on some ball $B_r$ to the above space of
$a$-harmonic functions that grow at most like a polynomial of degree $k$ has
the natural decay in the radius $r$ above some minimal radius $r_0$.
  Though motivated by stochastic homogenization, the contribution of this paper
is of purely deterministic nature: We work under the assumption that for the
given realization $a$ of the coefficient field, the couple $(\phi,\sigma)$ of
scalar and vector potentials of the harmonic coordinates, where $\phi$ is the
usual corrector, grows sublinearly in a mildly quantified way. We then
construct "$k$th-order correctors" and thereby the space of $a$-harmonic
functions that grow at most like a polynomial of degree $k$, establish the
above excess decay and then the corresponding Liouville principle.
