We study the approximation of harmonic functions by means of harmonic polynomials in
two-dimensional, bounded, star-shaped domains. Assuming that the functions possess
analytic extensions to a δ-neighbourhood of the domain, we prove
exponential convergence of the approximation error with respect to the degree of the
approximating harmonic polynomial. All the constants appearing in the bounds are explicit
and depend only on the shape-regularity of the domain and on δ. We apply
the obtained estimates to show exponential convergence with rate O(exp(-b√N)), N being the number of degrees of
freedom and b > 0, of a hp-dGFEM discretisation of
the Laplace equation based on piecewise harmonic polynomials. This result is an
improvement over the classical rate O(exp(-b3√N)), and is due to the use of harmonic polynomial
spaces, as opposed to complete polynomial spaces.