0$ such that \[ \big\|\frac{\Delta_h f}{\abs{h}}\big\|_{L^p(K)} \leq C \] for all $h \in \mathbb{R}^n$, $0<2\abs{h} <\operatorname{dist}(K,\partial \Omega) $. Then $f \in W^{1,p}(K)$. \end{enumerate} \end{theorem} Moreover, in the case $p=\infty$ we derive from \cite[Theorem 5, Section 4.2.3]{EG} the following result. In the proof we use $\alpha$-th difference quotients introduced in \eqref{dif2}. \begin{proposition}\label{wkp} Let $\Omega \subset \mathbb{R}^n$ be an open set, $f:\Omega \to \mathbb{R}$ and $k \in \mathbb{N}$, $k \ne 0$. \begin{enumerate} \item Suppose, that for all $K \Subset \Omega$ and all multi-indices $\alpha \in \mathbb{N}^n$ such that $|\alpha| \leq k $ there exists $C_{K,\alpha} >0$ such that for all $t=(t_1,\ldots,t_n) \in \mathbb{R}^n$ with $t_i \ne 0$ and $2 \abs{\alpha_1 t_1 + \ldots + \alpha_n t_n} < \operatorname{dist}(K,\partial \Omega)$ it holds \begin{equation}\label{e56} \big\|\frac{\Delta^\alpha_h f }{t_1^{\alpha_1} \cdots t_n^{\alpha_n} }\big\|_{L^\infty (K)} \leq C_{K,\alpha}, \end{equation} where $\Delta^\alpha_h$ is defined in \eqref{dif}, $h = et:=(t_1 e_1, t_2 e_2,\ldots,t_n e_n)$ and $e=(e_1,\ldots,e_n)$ is the standard basis of $\mathbb{R}^n$. Then $f \in W^{k,\infty}_{\rm loc}(\Omega)$. \item Suppose, that $f \in W^{k,\infty}(\Omega)$ and $K \Subset \Omega$. Then, for any $\alpha \in \mathbb{N}^n$, $|\alpha| \leq k $ there exists $C_{K,\alpha} >0$ such that for all $t=(t_1,\ldots,t_n) \in \mathbb{R}^n$ with $t_i \ne 0$ and $2 \abs{\alpha_1 t_1 + \ldots + \alpha_n t_n} < \operatorname{dist}(K,\partial \Omega)$ it holds \[ \big\|\frac{\Delta^\alpha_h f }{t_1^{\alpha_1} \cdots t_n^{\alpha_n}} \big\|_{L^\infty (K)} \leq C_{K,\alpha} \norm{f}_{W^{k,\infty}(\Omega)}. \] \end{enumerate} \end{proposition} \begin{proof} Fix $K \Subset \Omega$, $\alpha \in \mathbb{N}^n$, $|\alpha| \leq k$. We are going to show, that there exists an $\alpha$-th weak derivative of $f$ in $L^\infty(K)$. Let us fix a sequence $\{ t^k=(t_1^k, \ldots, t_n^k) \}_{k=1}^\infty$ such that $\left( t^k \right)^\alpha \to 0$. Observe, that by \eqref{e56} $\frac{\Delta^\alpha_h f }{t^\alpha}$ is bounded in $L^\infty(K)$, therefore by the Banach-Alaoglu Theorem, there exists a subsequence of $t^k$, still denoted in the same way, which has a weak-$*$ limit in $L^\infty(K)$ as $( t^k )^\alpha \to 0$. Denote this limit by $g^\alpha \in L^\infty(K)$. We need to show, that $g^\alpha$ is the $\alpha$-th weak derivative of $f$. Namely, we need to show that for any $\varphi \in C_c^\infty(K)$ the following holds: \begin{equation} \label{eq13} \int_{\mathbb{R}_n} \frac{\Delta^\alpha_h f (x)}{t^\alpha} \varphi(x) dx =(-1)^{|\alpha|} \int_{\mathbb{R}^n} f(x) \frac{\Delta^\alpha_{-h} \varphi (x)}{(-t)^\alpha} dx. \end{equation} This is an easy consequence of \begin{align*} \int_{\mathbb{R}^n} \frac{\Delta_h f(x) }{h} \varphi(x) dx &=\frac{1}{h} \int_{\mathbb{R}^n} (f(x+h)-f(x)) \varphi(x) dx \\ &= \frac{1}{h} \Big( \int_{\mathbb{R}^n} f(x) \varphi(x-h)dx - \int_{\mathbb{R}^n} f(x) \varphi(x) dx \Big) \\ & = \int_{\mathbb{R}^n} f(x) \frac{\varphi(x-h)-\varphi(x)}{h} dx \\ &=- \int_{\mathbb{R}^n} f(x) \frac{\Delta_{-h} \varphi(x)}{-h} dx. \end{align*} Notice, that for $t^\alpha \to 0$ the right-hand side of \eqref{eq13} converges to $\int_{\mathbb{R}^n} g^\alpha \varphi(x) dx$, and the left-hand side to $(-1)^{|\alpha|} \int_{\mathbb{R}^n} f(x) D^\alpha \varphi(x) dx$. Therefore $g^\alpha$ is the $\alpha$-th weak derivative of $f$, which completes the proof of the first assertion. The second assertion follows from $W^{k,\infty}(\Omega)= C^{k-1,1}(\Omega)$ and use of the Lipschitz condition. \end{proof} \section{Strongly harmonic functions on open subsets of $\mathbb{R}^n$} In this section we focus our attention on the class of strongly harmonic functions appearing in Definition~\ref{sh}. Let $(X,d,\mu)$ be a metric measure space with a Borel measure $\mu$. We denote by $\mathcal{H} (\Omega, d,\mu)$ the set of all strongly harmonic functions on an open domain $\Omega \subset X$. In what follows we will omit writing the set, metric or measure whenever they are clear from the context. The key results of this section are Proposition \ref{1p} and Theorem \ref{reg}. There, we show the Sobolev regularity for functions in $\mathcal{H} (\Omega,d,\mu)$ for the weighted Lebesgue measure $d \mu = w dx$ depending on the Sobolev regularity of weight $w$. The properties of strongly and weakly harmonic functions were broadly studied in \cite{H, J} and in~\cite{H1} in the setting of Carnot groups. Below, we list out some of those properties especially important for further considerations. \begin{proposition}[{\cite[Prop. 4.1]{H}}]\label{prop41} Suppose that $\mu$ is \emph{continuous with respect to metric}~$d$, i.e.\ for all $r>0$ and $x \in X$ it holds $ \lim_{d(x,y) \to 0} \mu\left( B(x,r)\triangle B(y,r) \right)=0$, where we denote by ${E \triangle F := (E \setminus F ) \cup (E \setminus F)}$ the symmetric difference of $E$ and $F$. Then $\mathcal{H}(\Omega,d,\mu)~\subset~C(\Omega)$. \end{proposition} Moreover, the Harnack inequality and the strong maximum principle hold for strongly harmonic functions as well as the local H\"older continuity and even local Lipschitz continuity under more involved assumptions, see \cite{H}. It is important to mention here that similar type of problems were studied for a more general, nonlinear mean value property by Manfredi-Parvainen-Rossi, Arroyo-Llorente, Ferrari-Liu-Manfredi and Ferrari-Pinamonti, see \cite{man,LL1,LL2,LL3,ferrari,FLM}. We know that $\mathcal{H}$ is a linear space, but verifying by using the definition whether some function satisfies the mean value property might be a complicated computational challenge. From that comes the need for finding a handy characterization of class $\mathcal{H}$, or some necessary and sufficient conditions for being strongly harmonic. Our goal is to characterize class $\mathcal{H}$ if $X = \mathbb{R}^n$ equipped with a distance $d$ induced by a norm and a weighted Lebesgue measure $d\mu=w dx$. From now on we a priori assume that a function $w \in L_{\rm loc}^1(\Omega)$ and $w>0$ almost everywhere in $\Omega$. Let us begin with noting that strongly harmonic functions in such setting are continuous. \begin{proposition} \label{prop4} Let $\Omega \subset \mathbb{R}^n$ be an open set. Then $\mathcal{H}(\Omega,d,w dx) \subset C(\Omega)$. \end{proposition} \begin{proof} Observe that $\mu(\partial B(x,r))= \int_{\partial B(x,r) } w(y) dy= 0$. Therefore, by \cite[Lemma 2.1]{J} measure $\mu$ is continuous with respect to metric. This completes the proof by Proposition \ref{prop41}. \end{proof} Let us observe that the proof of continuity of strongly harmonic functions works for all weights $w$. However, in order to show existence and integrability of weak derivatives we need to assume Sobolev regularity of $w$. \begin{proposition}\label{1p} Let $\Omega \subset \mathbb{R}^n$ be an open set, $d$ be a norm induced metric and a weight $w \in W^{1,p}_{\rm loc}(\Omega) \cap L^\infty_{\rm loc}(\Omega)$ for some $1

0$.
Then $ \mathcal{H}(\Omega,d,w ) \subset W^{l,\infty}_{\rm loc}(\Omega)$.
\end{theorem}
\begin{proof}
Let $u \in \mathcal{H}(\Omega,d,w)$ and $w$ be as in the assumptions. We will show that
$u \in W^{k,\infty}_{\rm loc}(\Omega)$ for every $k \leq l$ using the mathematical
induction with respect to $k$.
Let $k=1$ and $K,K',r,h,M$ be as in the proof of Proposition \ref{1p}.
The following is the consequence of \eqref{ew01}, \eqref{33} and \eqref{34},
\begin{align*}
\frac{\abs{\Delta_h u(x)}}{\abs{h}}
&\leq \norm{uw}_{_{L^\infty (K')}}
\Big( \frac{c_{n,d} r^{n-1}}{M} + \frac{\int_{B(x,r)} \abs{\Delta_h w}}{\abs{h} M^2} \Big) \\
&\leq \norm{uw}_{_{L^\infty (K')}} \Big(\frac{c_n r^{n-1}}{M}
+ \frac{C_{n,d} r^n }{M^2} \norm{ \frac{\Delta_h w}{\abs{h}}}_{L^\infty(K')} \Big)
\end{align*}
and is bounded by Proposition \ref{wkp}. Note, that this in particular means that
$u$ is Lipschitz.
Now let $k>1$ and assume that $u \in W^{k-1,\infty}_{\rm loc} (\Omega)$.
We consider the $\alpha$-th order difference quotient of $u$ for $|\alpha| =k$.
Let $t\in \mathbb{R}^n$ be such that ${|\alpha_1t_1 + \ldots + \alpha_n t_n|< \frac{r}{2k}}$
and define $h = (t_1 e_1, \ldots, t_n e_n)$. Formula \eqref{ee0}, Proposition \ref{pq}
and the related discussion applied to \sloppy ${f(x) = \int_{B(x,r)} uw}$ and
$g(x) = \int_{B(x,r)} w$ allows us to reduce the discussion to showing that
\[
\binom{\alpha}{\beta} \frac{ (-1)^m m!
\big( \Delta_h^{\alpha - \beta} \int_{B(x+\beta h,r)} uw \big)
\prod_{i=0}^m \big( \Delta_h^{\beta^i} \int_{B(x,r)} w \big)}
{t^\alpha \big( \int_{B(x,r)} w \big)^{m+1}}
\]
is bounded for any $\beta^1 + \ldots + \beta^m = \beta \leq \alpha$.
To show this we only need to show boundedness of the expression
\[
\frac{ \big( \Delta_h^{\alpha - \beta} \int_{B(x+\beta h,r)} uw \big)
\prod_{i=0}^m \big( \Delta_h^{\beta^i} \int_{B(x,r)} w \big)}{t^\alpha },
\]
since the term in the denominator is bounded from below by
\[
\Big( \int_{B(x,r)}w\Big)^{m+1} =\mu(B(x,r))^{m+1} \geq M^{m+1}
\]
and the rest of terms are constant. Let us observe that
\[
\frac{\Delta_h^{\beta^i} \int_{B(x,r)} w }{t^{\beta^i}}
= \int_{B(x,r)} \frac{\Delta_h^{\beta^i} w }{t^{\beta^i}}
\]
is bounded in $L^\infty(K)$ due to $w \in W^{1,\infty}_{\rm loc}(\Omega)$
and the second part of Proposition \ref{wkp}. Moreover, the upper bound is
constant multiple $\norm{w}_{W^{l,\infty}(K')}$.
Therefore, we only need to show boundedness of the term
$t^{\beta -\alpha}\Delta_h^{\alpha - \beta} \int_{B(x+\beta h,r)} uw $.
We only have to deal with the case $\beta = 0$. Indeed, observe, that for
$|\beta|>0$ the order of $(\alpha-\beta)$-th difference quotient is
$\abs{\alpha - \beta} \leq k-1$. Therefore, it is bounded due to the inductive
assumption that $uw \in W^{k-1,\infty}$. Observe that there exists $j$ such that
$\alpha_j \geq 1$ and so all components of vector $\alpha - e_j$ are non-negative natural.
Therefore the operator $\Delta_h^{\alpha-e_j}$ is well-defined,
$\Delta_h^\alpha =\Delta_{t_je_j} \circ \Delta_h^{\alpha - e_j}$, and hence
\begin{align*}
\frac{1}{\abs{h}^\alpha}\Big| \Delta^\alpha_h \int_{B(x,r)} uw\Big|
& = \frac{1}{\abs{h}^\alpha} \abs{ \Delta_{t_j e_j} \int_{B(x,r)}
\Delta^{\alpha - e_j}_h uw(y) dy} \\
& = \frac{1}{\abs{h}^\alpha} \Big|\int_{B(x+t_j e_j,r)}
\Delta^{\alpha-e_j}_h uw(y) dy -\int_{B(x,r)} \Delta^{\alpha-e_j}_h uw(y) dy \Big|\\
& \leq \frac{1}{\abs{t_j}} \int_{B(x+t_j e_j,r) \triangle B(x,r) }
\frac{\big|\Delta^{\alpha-e_j}_h uw(y) \big|}{\abs{h}^{\alpha-e_j}}dy\\
&\leq c_n r^{n-1}
\big\|\frac{ \Delta^{\alpha-e_j}_h uw }{\abs{h}^{\alpha-e_j}}\big\|_{L^\infty(K')},
\end{align*}
which is bounded by the regularity assumption on both $u$ and $w$ (in the last
estimate we have also used \eqref{35}). Therefore, we conclude that
$u\in W^{l,\infty}_{\rm loc}(\Omega)$, which completes the proof.
\end{proof}
\subsection{Historical background}
In what follows we are interested in extending results by Flatto \cite{F2,F3},
Friedman-Litmann \cite{F5}, Bose \cite{B2,B3,B4} and Zalcman \cite{F}.
Below, we briefly discuss these results. According to our best knowledge,
the investigation in this area originate from a work by Flatto \cite{F2}.
He considered functions with the following property:
Let us fix an open set $\Omega \subset \mathbb{R}^n$ and a bounded set
$K \subset \mathbb{R}^n$. Moreover, let $\mu$ be a probabilistic measure on $K$
such that all continuous functions on $K$ are $\mu$-measurable and for all
hyperplanes $V \subset \mathbb{R}^n$ it holds that $ \mu(K \cap V) <1$,
i.e.\ $\mu$ is not concentrated on a hyperplane. We will say that a continuous
function $u \in C(\Omega)$ has the mean value property in the sense of Flatto, if
\begin{equation}\label{fmean}
u(x) = \int_K u(x+ry) d\mu(y)
\end{equation}
for all $x \in \Omega$ and radii $r>0$ such that
$x +r \cdot K := \{x+r y: y \in K \} \subset \Omega$.
Let us observe that for $K = B(0,1)$ a unit ball in a given norm induced metric
$d$ and $\mu$ being the normalized Lebesgue measure on $K$, property \eqref{fmean}
is equivalent to the strong harmonicity of $u$ in $\Omega$ by the following formula
\begin{equation}\label{obs}
u(x)= -\hspace{-0.38cm}\int_{B(x,r)} u(z) dz
= -\hspace{-0.38cm}\int_{B(0,1)} u(x+r y ) dy
= \int_K u(x+r y ) d\mu(y).
\end{equation}
This holds exactly for homogeneous and translation invariant metrics, because only then
\[
B(x,r)=x+r \cdot B(0,1) =\{x+r y: y \in B(0,1) \}.
\]
For such distance functions one can obtain any ball $B(x,r)$ from $B(0,1)$
by using the change of variables $ y=\frac{z-x}{r}$. In relation to homogeneous
and translation invariant distance let us recall the following lemma, which
is likely a part of the mathematical folklore. However, in what follows we will
not appeal to this observation.
\begin{lemma} \label{lem2}
If $d$ is a translation invariant and homogeneous metric on $\mathbb{R}^n$,
then there exists a norm $\| \cdot \|$ on $\mathbb{R}^n$ such that for all
$x,y \in \Omega$ it holds that $d(x,y) = \| x-y \| $.
\end{lemma}
We present also a characterization of all such metrics on $\mathbb{R}^n$ using
the Minkowski functional, see~\cite{OO}. Recall, that a set
$K\subset \mathbb{R}^n$ is symmetric if $-y \in K$ for every $y \in K$.
For any nonempty convex set $K$ we consider the Minkowski functional.
\begin{lemma}[{\cite[p.54 ]{OO}}] \label{mink}
Suppose that $K$ is a symmetric convex bounded subset of $\mathbb{R}^n$, containing
the origin as an interior point. Then, its Minkowski functional $n_K$ defines
a norm on $\mathbb{R}^n$. Moreover, if $\|\cdot \|$ is a norm on $\mathbb{R}^n$,
then the Minkowski functional $n_K$, where $K$ is a unit ball with respect to
$\| \cdot \|$, is equal to that norm.
\end{lemma}
\begin{example}\label{uwmink} \rm
Among many examples of norm induced metrics on $\mathbb{R}^n$ are $l^p$ distances
for $1\leq p \leq \infty$. Moreover, let us fix numbers $a_i>0$ for $i=1,\ldots,n$,
set $a:=(a_1,\ldots,a_n)$ and $1\leq p< \infty$ and define
\[
\| x\|_p^a := \Big( \sum_{i=1}^n \Big( \frac{|x_i|}{a_i} \Big)^p \Big)^{1/p} .
\]
In case $p=2$ all balls with respect to $\| \cdot \|_p^a$ are ellipsoids with
the length of semi-axes equal to $a_i$ in $x_i$'s axes direction respectively.
Let us observe that by Lemma \ref{mink} there is the injective correspondence
between norms on $\mathbb{R}^n$ and a class of all symmetric convex open bounded
subsets $K$ of $\mathbb{R}^n$. More specifically, every $K$ defines a norm on
$\mathbb{R}^n$ through the Minkowski functional and vice versa, given a norm on
$\mathbb{R}^n$ the unit ball $B(0,1)$ is a symmetric convex open bounded set,
therefore provides an example of $K$. This can be expressed in one more way,
namely that all norms can be distinguished by their unit balls, so to construct
a norm we only need to say what is its unit ball. Therefore, further examples of
norms can be constructed for any $n$-dimensional symmetric convex polyhedron $K$.
All balls with respect to $n_K$ will be translated and dilated copies of $K$.
\end{example}
The formula \eqref{obs} is true only if the measure of a ball scales with the power
$n$ of its radius, the same which appears in the Jacobian from the change of
variables formula. This is true only for measures which are constant multiples of
the Lebesgue measure. Note that \eqref{fmean} does not coincide, in general,
with the mean value property presented in our work, since the Flatto's mean
value is calculated always with respect to the same fixed reference set $K$
and measure $\mu$, whose support is being shifted and scaled over $\Omega$.
Whereas, in Definition \ref{sh} the measure is defined on the whole space,
and as $x$ and $r$ vary, the mean value is calculated with respect to different
weighted measures. Indeed, from the point of view of Flatto, the condition from
Definition \ref{sh} reads
\[
u(x) = \int_{X} u(y) \frac{d \mu |_{B(x,r)}}{\mu(B(x,r))}.
\]
This mean value property cannot be written as an integral with respect to one fixed
measure for different pairs of $x$ and $r$, even when \eqref{obs} holds.
Flatto discovered that functions satisfying \eqref{fmean} are solutions to a
second order elliptic equation, see \cite{F2}. However, from the point of view
of our discussion, more relevant is the following later result.
\begin{theorem}[{Friedman-Littman \cite[Theorem 1]{F5}}]\label{F-L}
Suppose that $u$ has property \eqref{fmean} in $\Omega \subset \mathbb{R}^n$.
Then $u$ is analytic in $\Omega$ and satisfies the following system of partial
differential equations
\begin{equation}\label{system}
\sum_{|\alpha|=j} A_\alpha D^\alpha u =0 \quad \text{for } j=1,2,\ldots
\end{equation}
The coefficients are
$A_\alpha := \binom{|\alpha|}{\alpha} \int_K x^\alpha d\mu(x)$ and
are moments of measure $\mu$.
Moreover, any function $u \in C^\infty (\Omega)$ solving system \eqref{system}
is analytic and has property \eqref{fmean}.
\end{theorem}
Theorem gives full characterization of $\mathcal{H}(\Omega,d)$ for $d$ being induced by a norm.
\cite[Theorem 3.1]{F2} states that all functions having property \eqref{fmean}
are harmonic with respect to variables obtained from $x$ by using an orthogonal
transformation and dilations along the axes of the coordinate system.
On the other hand the proof of Theorem \ref{F-L} shows that the equation in
system \eqref{system} corresponding to $j=2$ is always elliptic with constant
coefficients from which the analyticity follows.
Flatto as well as Friedman and Littman described in their works the space of
functions possessing property \eqref{fmean}. We present appropriate results below.
\begin{proposition}[{Friedman-Littman \cite[Theorem 2]{F5}}]\label{FL2}
The space of solutions to system \eqref{system} is finitely dimensional if and only
if the system of algebraic equations $\sum_{|\alpha|=j} A_\alpha z^\alpha = 0$
for $j=1,2,\ldots$ has the unique solution $z=(z_1,\ldots,z_n)=0$, where
$z_i \in \mathbb{C}$.
\end{proposition}
\begin{remark}\label{uw6} \rm
From the proof of Proposition \ref{FL2} it follows that if there exists
a nonpolynomial solution to~\eqref{system}, then the solution space is infinitely
dimensional. If the dimension is finite, then all strongly harmonic functions are polynomials.
\end{remark}
A different approach to the mean value property and its consequences was
studied by Bose, see \cite{B2,B3,B4}. He considered strongly harmonic functions
on $\Omega \subset \mathbb{R}^n$ equipped with non-negatively weighted measure
$\mu = w dx$, for a weight $w \in L_{\rm loc}^1(\Omega)$ being a.e. positive
in $\Omega$ and only a metric $d$ induced by the $l_2$-norm.
Under the higher regularity assumption of weight $w$, Bose proved the following result.
\begin{theorem}[{Bose \cite[Thm.\ 1]{B4}}]\label{twbo}
If there exists $l \in \mathbb{N}$ such that $w \in C^{2l+1}(\Omega)$, then every
$u \in \mathcal{H}(\Omega,w)$ solves the following system of partial differential equations
\begin{equation}\label{system2}
\Delta u \Delta^j w + 2 \nabla u \nabla \left( \Delta^j w \right)=0, \quad
\text{for } j=0,1,\ldots,l,
\end{equation}
where $\Delta^j$ stands for the $j$th composition of the Laplace operator $\Delta$.
If $w$ is smooth, then equations \eqref{system2} hold true for all $j \in \mathbb{N}$.
\end{theorem}
The converse is not true for smooth weights in general, see counterexamples in
\cite[p. 479]{B2}. Furthermore, Bose proved in \cite{B4} the following result,
by imposing further assumptions on $w$.
\begin{proposition}[Bose]\label{bosod}
Let $l \in \mathbb{N}$ and $w \in C^{2l}(\Omega)$. Suppose that there exist
$a_0,\ldots,a_{l-1} \in \mathbb{R}$ such that
\[
\Delta^l w = a_0 w + a_1 \Delta w + \ldots + a_{l-1} \Delta^{l-1} w.
\]
Then any solution $u$ to \eqref{system2} for $j=0,1,\ldots,l-1$ is strongly harmonic,
that is $u \in \mathcal{H}(\Omega,w)$.
\end{proposition}
The following result by Bose contributes to the studies of the dimension of the space
$\mathcal{H}(\Omega,l^2,w)$ under certain additional assumption on the weight
(in particular, assuming that $w$ is an eigenfunction for the laplacian).
\begin{proposition}[{Bose \cite[Corollary 2]{B2}}]
Suppose that $\Omega \subset \mathbb{R}^n$ for $n>1$, $w \in C^2(\Omega)$
and there exists $\lambda \in \mathbb{R}$ such that $ \Delta w = \lambda w$.
Then $\operatorname{dim} \mathcal{H}(\Omega,w) = \infty$.
\end{proposition}
\begin{remark}\label{rem1} \rm
Before we present the proof of Theorem \ref{twsys} let us discuss the equations
of system \eqref{system3}. First of all, by Remark \ref{uwmink} we know that the
set $B(0,1)$ is symmetric with respect to the origin.
If $|\alpha|$ is an odd number, then $x^\alpha$ is an odd function, hence $A_\alpha = 0$.
Therefore only evenly indexed equations of \eqref{system3} are nontrivial, although
we will prove them for all $j \leq 2l$. In fact, the proof of Theorem \ref{twsys}
can be applied to functions with the mean value property over any compact set
$K \subset \mathbb{R}^n$, which does not necessarily need to be a unit ball with
respect to a norm on $\mathbb{R}^n$, i.e. to functions with the following property
\[
u(x) = \frac{1}{\int_K w(x+ry) dy}\int_{K} u(x+ry)w(x+ry) dy,
\]
which holds for all $x \in \Omega$ and radii $0