Cubic harmonics and Bernoulli numbers

The functions satisfying the mean value property for an n-dimensional cube are determined explicitly. This problem is related to invariant theory for a finite reflection group, especially to a system of invariant differential equations. Solving this problem is reduced to showing that a certain set of invariant polynomials forms an invariant basis. After establishing a certain summation formula over Young diagrams, the latter problem is settled by considering a recursion formula involving Bernoulli numbers.