\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{amssymb} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2020 (2020), No. 08, pp. 1--26.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2020 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2020/08\hfil Characterization of mean value harmonic functions] {Characterization of mean value harmonic functions on norm induced metric measure spaces with weighted Lebesgue measure} \author[A. Kijowski \hfil EJDE-2020/08\hfilneg] {Antoni Kijowski} \address{Antoni Kijowski \newline Institute of Mathematics, Polish Academy of Sciences, ul. \'Sniadeckich 8, 00-656 Warsaw, Poland} \email{akijowski@impan.pl} \thanks{Submitted March 21, 2019. Published January 14, 2020.} \subjclass{31C05, 35J99, 30L99} \keywords{Harmonic function; mean value property; metric measure space; \hfill\break\indent Minkowski functional; norm induced metric; Pizzetti formula; weighted Lebesgue measure} \begin{abstract} We study the mean-value harmonic functions on open subsets of~$\mathbb{R}^n$ equipped with weighted Lebesgue measures and norm induced metrics. Our main result is a necessary condition stating that all such functions solve a certain homogeneous system of elliptic PDEs. Moreover, a converse result is established in case of analytic weights. Assuming the Sobolev regularity of the weight $w \in W^{l,\infty}$ we show that strongly harmonic functions are also in~$W^{l,\infty}$ and that they are analytic, whenever the weight is analytic. The analysis is illustrated by finding all mean-value harmonic functions in~$\mathbb{R}^2$ for the $l^p$-distance ${1 \leq p \leq \infty}$. The essential outcome is a certain discontinuity with respect to $p$, i.e.\ that for all $p \ne 2$ there are only finitely many linearly independent mean-value harmonic functions, while for $p=2$ there are infinitely many of them. We conclude with the remarkable observation that strongly harmonic functions in~$\mathbb{R}^n$ possess the mean value property with respect to infinitely many weight functions obtained from a given weight. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{definition}[theorem]{Definition} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{remark}[theorem]{Remark} \newtheorem{example}[theorem]{Example} \allowdisplaybreaks \newcommand{\norm}{\|#1\|} \newcommand{\abs}{|#1|} \section{Introduction} Analysis on metric spaces has been intensively developed through the previous two decades. Studies of such researchers as Cheeger, Haj{\l}asz, Heinonen, Koskela and Shanmugalingam brought new light to a notion of the gradient in metric measure spaces. One of many important notions of this area is a counterpart of a harmonic function on metric measure spaces being a minimizer of the Dirichlet energy. Recently, there has been a new approach to this topic by using the mean value property. Such an approach is much easier to formulate, than the variational one, because it does not require the notion of the Sobolev spaces on metric measure spaces. Strongly and weakly harmonic were introduced in \cite{H,J} by Adamowicz, Gaczkowski and G\'orka. Authors developed the theory of such functions providing e.g.\ the Harnack inequality, the H\"older and Lipschitz regularity results and studying the Perron method. Nevertheless, many questions remain unanswered, including the one on the relation between minimizers of the Dirichlet energy and mean value harmonic functions. In order to understand this class of functions in the abstract metric setting one needs to investigate their properties in the classical setting of Euclidean domains, or in the wider class of Riemannian manifolds. Recall, that by a \emph{metric measure space} we mean metric space $(X,d)$ equipped with Borel regular measure $\mu$, which assigns to every ball a positive and finite value. In this setting we introduce the following class of functions. \begin{definition}[{\cite[Definition 3.1]{H}}] \label{sh} \rm Let $\Omega \subset X$ be an open set. We say that a locally integrable function $u: \Omega \to \mathbb{R}$ is \emph{strongly harmonic} in $\Omega$ if for all open balls $B(x,r) \Subset \Omega$ it holds $u(x) = -\hspace{-0.38cm}\int_{B(x,r)} u(y) d\mu(y) := \frac{1}{\mu(B(x,r))} \int_{B(x,r)} u(y) d\mu(y).$ We call a radius $r>0$ \emph{admissible} at some $x \in \Omega$ whenever $B(x,r) \Subset \Omega$. The space of all strongly harmonic functions in $\Omega$ is denoted by $\mathcal{H}(\Omega,d,\mu)$. We omit in this notation writing the set, metric or measure whenever they are clear from the context. \end{definition} The main object of this work is a characterization of strongly harmonic functions on a certain class of metric measure spaces. Namely, we consider an open subset $\Omega \subset \mathbb{R}^n$ equipped with a weighted Lebesgue measure $d\mu = w dx$, $w \in L^1_{\rm loc}(\Omega)$, $w>0$ a.e.\ and a \emph{norm induced metric} $d$, i.e.\ $d$ is a metric on $\Omega$ such that there exists a norm $\norm{\cdot}:\mathbb{R}^n \to [0,\infty)$ and for every $x,y \in \Omega$ it holds $d(x,y)=\norm{x-y}$. Bose, Flatto, Friedman, Littman, Zalcman studied the mean value property in the Euclidean setting, see \cite{B2,B3,B4,F2,F3,F1,F5,F}. We extended their results with our main result, see Theorem \ref{twsys} below. It generalizes results in \cite{F5} (see Theorem \ref{F-L} below) and in \cite{B4} (see Theorem \ref{twbo} below) in the following ways: \begin{enumerate} \item[(1)] we consider general metric functions induced by a norm, not necessarily the Euclidean one, \item[(2)] we allow more general measures, i.e. the weighted Lebesgue measures $d\mu = w dx$, \end{enumerate} Throughout this article we use the multi-index notation: $\alpha = (\alpha_1,\ldots,\alpha_n) \in \mathbb{N}^n$, $|\alpha| = \alpha_1 + \ldots + \alpha_n$. For more information see Appendix A in the Evans' book \cite{E2}. \begin{theorem}\label{twsys} Let $\Omega \subset \mathbb{R}^n$ be an open set. Let further $(\Omega,d,\mu)$ be a metric measure space equipped with a norm induced metric $d$ and a weighted Lebesgue measure $d\mu=wdx$, $w \in L^1_{\rm loc}(\Omega)$, $w>0$ a.e.\ Suppose that the weight $w \in W^{2l,\infty}_{\rm loc}(\Omega)$ for some given $l \in \mathbb{N}$, $l>0$. Then every function $u \in \mathcal{H}(\Omega,d,w dx)$ is a weak solution to the following system of partial differential equations \begin{equation}\label{system3} \sum_{|\alpha|=j} A_\alpha \left( D^\alpha (uw) - uD^\alpha w \right) =0, \quad \text{for }j=2,4,\ldots,2l. \end{equation} The coefficients $A_\alpha$ are defined as $A_\alpha:= \binom{|\alpha|}{\alpha} \int_{B(0,1)} x^\alpha dx = \frac{|\alpha|!}{\alpha_1! \cdots \alpha_n!} \int_{B(0,1)} x_1^{\alpha_1}\cdots x_n^{\alpha_n} dx,$ where $B(0,1)$ is a unit ball in metric $d$. \end{theorem} Let us briefly compare Theorem \ref{twsys} with $d=l^2$ to Bose's results \cite{B2,B3,B4}. In order to prove the necessary condition Theorem \ref{twbo} for being strongly harmonic Bose assumes the regularity of weight $w \in C^{2l-1}(\Omega)$, whereas our methods for showing Theorem \ref{twsys} require that $w \in W^{2l,\infty}_{\rm loc}(\Omega) = C^{2l-1,1}(\Omega)$. Nevertheless, if $d=l^2$ we retrieve the same system of PDEs as Bose, however this observation needs additional calculations presented in Section 4.1. On the other hand, to prove the sufficient condition for being strongly harmonic Bose assumes that the weight $w$ is an generalized eigenfunction of the laplacian, see Proposition~\ref{bosod}. In Theorem \ref{twod} we assume analyticity of weight $w$ to prove the sufficient condition. Our assumption is more general than Bose's, which is illustrated by Lemma \ref{lembo}. To prove Theorem \ref{twsys} we need to establish regularity results which are stated as Proposition~\ref{1p} and Theorem \ref{reg}. Roughly speaking, Proposition \ref{1p} shows that if weight $w$ is locally bounded and in the space $W^{1,p}_{\rm loc}$, then all strongly harmonic functions are in $W^{1,p}_{\rm loc}$, while Theorem \ref{reg} says that if $w$ is in $W^{l,\infty}_{\rm loc}$, then all strongly harmonic functions are in $W^{l,\infty}_{\rm loc}$. The discussion demonstrating the way how Theorem \ref{twsys} generalizes Theorem \ref{twbo} requires computations. We present them after the proof of Theorem \ref{twsys}, in Section~4.1. Our second main result is the following converse to Theorem \ref{twsys}. \begin{theorem}\label{twod} Let $\Omega \subset \mathbb{R}^n$ be an open set and $(\Omega,d,\mu)$ be a metric measure space equipped with a norm induced metric $d$ and a weighted Lebesgue measure $d\mu = w dx$. Suppose that weight $w$ is analytic and positive in $\Omega$. Then, any solution $u$ to system of equations \eqref{system3} is strongly harmonic in $\Omega$. \end{theorem} Another, perhaps most surprising results are presented in Section 4 where we illustrate Theorem \ref{twsys} with the following observations: \emph{If $p\ne 2$ and $n=2$, then the space $\mathcal{H}(\Omega,l^p,dx)$ is spanned by 8 linearly independent harmonic polynomials.} We already know that for any $n\geq 1$ the space $\mathcal{H}(\Omega,l^2,dx)$ consists of all harmonic functions in $\Omega$, and is infinitely dimensional. The result describing $\dim \mathcal{H}(\Omega,l^p,dx )$ for $p \ne 2$ in dimension $n=3$ is due to \L ysik \cite{GL}, who computed it to be equal to 48. The problem for $n > 3$ is open. It is also worthy mentioning here, that the dimensions 8 for $n=2$ and 48 for $n=3$ coincide with the number of linear isometries of the normed space $(\mathbb{R}^n, l^p)$, which is $2^n n!$ and is computed in \cite{LI}. For more information see Sections 4.2, 4.3 and 5.1. \subsection*{Organization of this article} In the preliminaries we introduce basic notions and definitions, which will be essential in further parts of the paper. The difference quotients characterization of Sobolev spaces is recalled and the formula for difference quotients of a quotient of two functions is developed. In Section 3 we present a historical sketch of the studies of the mean value property ending with the proof of Theorem \ref{twsys}. Moreover, by assuming the Sobolev regularity of weights, we prove in Theorem \ref{reg} that strongly harmonic functions are in the Sobolev space of the same order as the weight, see also Proposition~\ref{1p}. Further on we recall results of Flatto and Friedman-Littman concerning functions with the mean value property in the sense of Flatto (see \eqref{fmean} below) and compare them to strongly harmonic functions. Then, we recall a result of Friedman--Littman \cite{F5} which characterizes functions with the mean value property in the sense of Flato for the Lebesgue measure, but a metric not necessarily the Euclidean one. In fact we extend their proof to describe such functions via a system of PDEs. On the other hand, we present another approach studied by Bose \cite{B2,B3,B4}. He considered a mean value property on Euclidean balls for a weighted Lebesgue measure. We generalize both approaches in Theorem \ref{twsys} to the case of a weighted Lebesgue measure and a norm induced metric. We show that this case is the only one in which strongly harmonic functions coincide with those having mean value property in the sense of Flatto. In Section 4 we focus on $l^p$ metrics for $1\leq p \leq \infty$. Equations of system~\eqref{system3} are calculated explicitly with their coefficients $A_\alpha$. We show that there appear only two distinct cases: either $p=2$ and $\mathcal{H}(\Omega,l^2,dx)$ consists of all functions which Laplacian vanishes in $\Omega$, or $p\ne 2$ and there are only finitely many linearly independent strongly harmonic functions in the space $\mathcal{H}(\Omega,l^p,dx)$. Similar observations can be obtained in higher dimensions using our techniques. The last Section is devoted to proving Theorem \ref{twod}, a converse to Theorem \ref{twsys}. In order to complete that goal we recall the notion of generalized Pizzetti formula following Zalcman \cite{F}. Moreover, in Lemma \ref{lem4} we prove that equation for $j=2$ of \eqref{system3} is of the elliptic type. We use this fact to prove regularity of strongly harmonic functions, i.e. that all strongly harmonic functions are analytic whenever weight is analytic. We conclude Section 5 with applying Theorem \ref{twod} to obtain the following peculiar observation. Suppose that $u$ is strongly harmonic, weight $w$ is smooth and metric is Euclidean. Then, $u$ is strongly harmonic with respect to infinitely many weights obtained as compositions of the Laplacian on $w$, {i.e. $\Delta^l w$ for $l \in \mathbb{N}$.} \section{Preliminaries} In this section we introduce basic notions used in further parts of the work. Let $\mu$ be a measure on $\mathbb{R}^n$, set $A \subset \mathbb{R}^n$ be of positive measure $\mu(A)\, >0$, and $f: \mathbb{R}^n\to \mathbb{R}$ a measurable function. Then, the \emph{mean value of $f$} over set $A$ will be denoted by \begin{equation*} -\hspace{-0.38cm}\int_{A} f(x) d\mu(x) := \frac{1}{\mu(A)}\int_{A} f(x) d\mu(x). \end{equation*} We say that a function $u \in C^2(\Omega)$ is \emph{harmonic} in an open set $\Omega \subset \mathbb{R}^n$, if $\Delta u =0$ in $\Omega$. One of several properties of harmonic functions is the Gauss theorem stating that if $u$ is harmonic, then it has the mean value property with respect to the Lebesgue measure on all balls and spheres. There is an elegant converse relation between the mean value property and harmonicity brought by Hansen-Nadirashvili in \cite{A}: Let $\Omega$ be an open bounded subset of $\mathbb{R}^n$, $u \in C( \Omega) \cap L^\infty(\Omega)$ be such that for every $x \in \Omega$ there exists $0< r^x \leq \operatorname{dist}(x,\partial \Omega)$ with the property $u(x) = -\hspace{-0.32cm}\int_{B(x,r^x)} u(y) dy$. Then $u$ is harmonic in $\Omega$. The aforementioned relation between harmonicity and the mean value property leads to formulating a relaxed version of the strong harmonicity (cf. Definition \ref{sh}): Let $\Omega \subset X$ be an open set. We call any locally integrable function $u: \Omega \to \mathbb{R}$ \emph{weakly harmonic} in $\Omega$ if for all points $x \in \Omega$ there exists at least one radius \break $0 < r^x < \operatorname{dist}(x, \partial \Omega)$ with the following property $u(x) = -\hspace{-0.32cm}\int_{B(x,r^x)} u(y) d\mu(y)$. For further information about properties of weakly and strongly harmonic functions we refer to \cite{H,J}. Let us consider a function $f: \mathbb{R}^n \to \mathbb{R}$. For $x,h \in \mathbb{R}^n$ we define the \emph{difference} of $f$ at $x$ as follows \begin{equation}\label{d} \Delta_h f(x) := f(x+h) -f(x). \end{equation} Observe, that for any $h,h' \in \mathbb{R}^n$ difference operators $\Delta_h$ and $\Delta_{h'}$ commute and that $\Delta_0 f \equiv 0$. \begin{lemma}\label{lems} For a smooth function $f:\mathbb{R} \to \mathbb{R}$, its $k$-th derivative $f^{(k)}$ and $h_1,\ldots ,h_k \in \mathbb{R} \setminus \{0\}$ it holds $\frac{\Delta_{h_1} \circ \ldots \circ \Delta_{h_k} f(x)}{h_1 \cdots h_k} = f^{(k)}(x) + O( |h_1| + \ldots + |h_k|),$ as $|h_1| + \ldots + |h_k| \to 0$. \end{lemma} \begin{proof} We will show the assertion using the mathematical induction with respect to $k$. For $k=1$ and $h_1 \in \mathbb{R} \setminus \{ 0\}$ we apply the Taylor expansion theorem to obtain that $f(x+h_1)-f(x) = h_1 f'(x) + \frac{h_1^2}{2}f''(x) + o(h_1^2)$ as $h_1 \to 0$. Therefore $\frac{\Delta_{h_1} f(x)}{h_1} - f'(x) = \frac{h_1}{2}f''(x) + \frac{o(h_1^2)}{h_1} = O(|h_1|).$ For given $k>1$ and $h_1,\ldots,h_k \in \mathbb{R} \setminus \{ 0\}$ the inductive assumption reads \begin{align*} \frac{\Delta_{h_1} \circ \ldots \circ \Delta_{h_{k-1}} f(x)}{h_1 \cdots h_{k-1}} &= f^{(k-1)}(x) + O(|h_1| + \ldots + |h_{k-1}|). \end{align*} Observe, that the left-hand side is a smooth function with respect to $x$ and that derivatives commute with difference quotients in the following sense $\Big( \frac{\Delta_{h_1} \circ \ldots \circ \Delta_{h_{k-1}} f(x)}{h_1 \cdots h_{k-1}} \Big)' = \frac{\Delta_{h_1} \circ \ldots \circ \Delta_{h_{k-1}} f'(x)}{h_1 \cdots h_{k-1}}.$ Therefore, by the result for $k=1$ we have that \begin{align*} \frac{\Delta_{h_1} \circ \ldots \circ \Delta_{h_k} f(x)}{h_1 \cdots h_k} & = \frac{1}{h_k} \Delta_{h_k} \Big( \frac{\Delta_{h_1} \circ \ldots \circ \Delta_{h_{k-1}} f(x)}{h_1 \cdots h_{k-1}} \Big)\\ & = \frac{\Delta_{h_1} \circ \ldots \circ \Delta_{h_{k-1}} f'(x)}{h_1 \cdots h_{k-1}} + O(|h_k|). \end{align*} Applying the inductive assumption to $f'(x)$ we obtain that \begin{align*} \frac{\Delta_{h_1} \circ \ldots \circ \Delta_{h_{k-1}} f'(x)}{h_1 \cdots h_{k-1}} + O(|h_k|) &= f^{(k)}(x) + O(|h_1| + \ldots + |h_{k-1}|) +O(|h_k|)\\ & = f^{(k)}(x) + O( |h_1| + \ldots + |h_k|). \end{align*} From this, the assertion of the lemma follows. \end{proof} For $h=(h_1,\ldots ,h_n )$, $h_i \in \mathbb{R}^n$ and $\alpha=(\alpha_1,\ldots ,\alpha_n) \in \mathbb{N}^n$ the \emph{$\alpha$-th difference} of $f$ is defined as follows \begin{equation}\label{dif} \Delta_h^\alpha f(x) := (\Delta_{h_1})^{\alpha_1} \circ (\Delta_{h_2})^{\alpha_2} \circ \ldots \circ (\Delta_{h_n})^{\alpha_n} f(x) . \end{equation} We denote by the \emph{$\alpha$-th difference quotient} of $f$ the expression \begin{equation}\label{dif2} \frac{\Delta_h^\alpha f(x)}{h^\alpha} : =\frac{\Delta_h^\alpha f(x)}{h_1^{\alpha_1} \cdots h_n^{\alpha_n}}, \end{equation} whenever $h_i \ne 0$ for $i=1,\ldots,n$ such that $\alpha_i\ne 0$ (here we interpret the symbol $0^0$ to be equal to $1$). Formulas describing the difference quotients of a multiple and a quotient of two functions are similar to formulas describing their derivatives. Let us consider two functions $f,g: \mathbb{R}^n \to \mathbb{R}$ with $g>0$. In what follows we will need a representation of the $\alpha$-th difference quotient of $f/g$ in terms of difference operators applied to the nominator $f$ and the denominator $g$. Observe, that for $\alpha \in \mathbb{N}^n$ it holds a differences analogue of the Leibniz formula \begin{equation}\label{ee0} \Delta^\alpha_h \Big( \frac{f}{g} \Big) (x) = \sum_{\beta \leq \alpha} \binom{\alpha}{\beta} \Delta_h^{\alpha - \beta} f(x + \beta h) \Delta_h^\beta \Big( \frac{1}{g} \Big)(x), \end{equation} where $\beta h := (\beta_1 h_1, \ldots ,\beta_n h_n)$, notation $\beta \leq \alpha$ means $\beta_i \leq \alpha_i$ for all $i=1,\ldots,n$ and we denote by $\binom{\alpha}{\beta}=\binom{\alpha_1}{\beta_1} \cdots \binom{\alpha_n}{\beta_n}$. Notice, that when using \eqref{ee0} we only need to express the $\beta$-th difference quotient of the function $1/g$ in terms of difference quotients of $g$. To reach that goal we use a discrete variant of the Fa\'a di Bruno formula developed in \cite{FAA}, from which one can derive the following result. \begin{proposition}\label{pq} Let $\beta \in \mathbb{N}^n$, $x \in \mathbb{R}^n$, $h =(h_1,\ldots,h_n )$, $h_i \in \mathbb{R}^n \setminus \{0\}$ and $g:\mathbb{R}^n \to \mathbb{R}$ be a positive continuous function. Then \begin{equation}\label{32} \Delta^\beta_h \Big(\frac{1}{g}\Big) (x) = \sum_{ \beta^1 + \ldots + \beta^m = \beta} \frac{(-1)^m m!}{g(x)^{m+1}} \Delta_h^{\beta^1} g(x) \cdots \Delta_h^{\beta^m} g(x) + Err(x,h), \end{equation} where we sum with respect to $\beta^i \in \mathbb{N}^n \setminus \{0\}$ for $i=1,\ldots,m$, $m\in \mathbb{N}$, $m>0$ for $|h_1|+ \ldots + |h_n|$ small enough. The expression $Err(x,h)$ contains terms of order higher than $|\beta|$ in the sense defined in \cite{FAA} and is precisely described in \cite[Theorem 1.4]{FAA}. \end{proposition} Before we present the proof of Proposition \ref{pq} we want to give the reader some intuition by proving Proposition \ref{pq} for $n=1$ and $|\beta|=1$. Consider $h_1=:h \in \mathbb{R} \setminus \{ 0 \}$ and write \begin{align*} \Delta_{h_1}^\beta \Big(\frac{1}{g}\Big) (x) &= \Delta_h \Big(\frac{1}{g}\Big) (x) = \frac{1}{g(x+h)} - \frac{1}{g(x)} \\ &= \frac{\frac{1}{g(x+h)} - \frac{1}{g(x)}}{g(x+h)-g(x)} ( g(x+h) - g(x)) \\ &= \Big( \frac{d}{dy} \frac{1}{y} \big|_{y=g(x)} + O(\abs{\Delta_h(g(x))}) \Big) \Delta_h g(x) \\ &= \frac{-1}{g(x)^2} \Delta_h g(x) + O(\abs{\Delta_h(g(x))})\Delta_h g(x), \end{align*} where passing from the first to the second line we used Lemma \ref{lems} with $f(y)=1/y$. \begin{proof}[Proof of Proposition \ref{pq}] Recall that by \eqref{d}, $$\Delta_{\Delta^{\beta^i}_h g(x)} \frac{1}{g(x)} = \frac{1}{g\big( x+\Delta^{\beta^i}_h g(x) \big)} - \frac{1}{g(x)}.$$ We apply \cite[Theorem 1.3]{FAA} to the composition of function $f(x)=\frac1x$ with $g: \mathbb{R}^n \to \mathbb{R}$, $\Delta_h^\beta \Big(\frac{1}{g}\Big)(x) = \sum_{\beta^1 + \ldots + \beta^m = \beta} \Delta_{\Delta_h^{\beta^1} g(x)} \circ \ldots \circ \Delta_{\Delta_h^{\beta^m} g(x)} \frac{1}{g(x)} + Err'(x,h),$ where $Err'(x,h)$ stands for the error term in \cite[Theorem 1.3]{FAA}. Let us calculate the following \begin{align*} &\sum_{ \beta^1 + \ldots + \beta^m = \beta} \Delta_{\Delta_h^{\beta^1} g(x)} \circ \ldots \circ \Delta_{\Delta_h^{\beta^m} g(x)} \frac{1}{g(x)} \\ &= \sum_{ \beta^1 + \ldots + \beta^m = \beta} \frac{\Delta_{\Delta_h^{\beta^1} g(x)} \circ \ldots \circ \Delta_{\Delta_h^{\beta^m} g(x)} \big( \frac{1}{g(x)} \big)}{\Delta_h^{\beta^1} g(x) \cdots \Delta_h^{\beta^m} g(x) } \cdot \Delta_h^{\beta^1} g(x) \cdots \Delta_h^{\beta^m} g(x), \end{align*} where we omit writing these terms in the sum, for which at least one of $\Delta_h^{\beta^i} g(x)=0$ for $i=1,\ldots,m$, since then the corresponding term vanishes. Now let us observe, that the quotient appearing under the sum is in fact the $m$-th order difference quotient of the function $\frac1x$ at point $g(x)$ with increments $\Delta^{\beta^i}_h g(x)$, which all by the continuity assumption of $g$ tend to $0$ as $h^{\beta^i} \to 0$. Therefore, by Lemma \ref{lems} we obtain \begin{align*} \Delta_h^\beta \Big(\frac{1}{g}\Big)(x) &= \sum_{ \beta^1 + \ldots + \beta^m = \beta} \Big( \frac{d^m}{dy^m} \frac{1}{y} \big|_{y=g(x)} + O \Big( |\Delta_h^{\beta^1} g(x)| + \ldots + |\Delta_h^{\beta^m} g(x)| \Big) \Big) \\ &\quad\times \Delta_h^{\beta^1} g(x) \cdots \Delta_h^{\beta^m} g(x) + Err'(x,h) \\ &= \sum_{ \beta^1 + \ldots + \beta^m = \beta} \frac{(-1)^m m!}{g(x)^{m+1}} \Delta_h^{\beta^1} g(x) \cdots \Delta_h^{\beta^m} g(x)\\ &\quad + \sum_{ \beta^1 + \ldots + \beta^m = \beta} O \left( \abs{\Delta_h^{\beta^1} g(x)} + \ldots + \abs{\Delta_h^{\beta^m} g(x)} \right) \Delta_h^{\beta^1} g(x) \cdots \Delta_h^{\beta^m} g(x) \\ &\quad + Err'(x,h) \\ & =\sum_{ \beta^1 + \ldots + \beta^m = \beta} \frac{(-1)^m m!}{g(x)^{m+1}} \Delta_h^{\beta^1} g(x) \ldots \Delta_h^{\beta^m} g(x) + Err(x,h), \end{align*} which completes the proof. \end{proof} An outcome of the above discussion is that we can represent the $\alpha$-th difference quotient of $f/g$ as the sum of fractions whose numerators, apart from constants, consist only of terms $\Delta^{\beta-\alpha}_h f(x+ \beta h)$, $\Delta_h^{\beta^j} g(x)$ and their products for some $\beta^1 + \ldots + \beta^m = \beta \leq \alpha$. Furthermore, the operator $\Delta_h$ appears in each of these numerators exactly $|\alpha|$-times, which can be justified by calculating the sum $\sum_{j=1}^m \abs{\beta_j} + \abs{\alpha-\beta}= |\alpha|$. We use difference quotients to prove regularity of strongly harmonic functions in Theorem \ref{reg}. Therefore, we gather below a characterization of Sobolev functions via difference quotients, recall \eqref{d}. \begin{theorem}[{\cite[Theorem 3, p. 277]{E2}}] \label{w1p} Let $\Omega \subset \mathbb{R}^n$ be an open set. \begin{enumerate} \item Suppose that $1 \leq p < \infty$, $f \in W^{1,p}(\Omega)$. Then for each $K \Subset \Omega$ $\big\|\frac{\Delta_{h} f}{\abs{h}}\big\|_{L^p(K)} \leq C \norm{\nabla f}_{L^p(\Omega)},$ for some constant $C>0$ and all $h \in \mathbb{R}^n$, $0<2\abs{h} <\operatorname{dist}(K,\partial \Omega)$. \item Suppose that $10$ 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 $10$. 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 $00$. Furthermore, let $B'$ be a ball concentric with $B$ with $\varepsilon$ distance from $\partial \Omega$. We redefine $u$ and $w$ in the following way \begin{equation*} \bar{u}(x)= u(x) \chi_{B'}(x) \quad \bar{w}(x) = w(x) \chi_{B'}(x). \end{equation*} The function $\bar{u}$ and the weight $\bar{w}$ are both in $W^{2l,\infty}(B)$ since $B \Subset B'$. Let $\varphi \in C^\infty_0 (B)$. Then for all $x \in B$, $y \in B(0,1)$ and $02$ since the equation for $j=2$ is described in \eqref{eq11}. We examine the differential operator $R_j:=\sum_{|\alpha|=4} A_\alpha D^\alpha$. We already showed that for $p=2$ operator $R_2$ is equal to $\Delta$ up to a multiplicative constant. Recall formula \eqref{ani}: $A_\alpha = \binom{|\alpha|}{\alpha} \Big( \frac{2}{p} \Big)^2 \frac{\prod_{i=1}^2 \Gamma\big( \frac{\alpha_i +1}{p} \big)} {\Gamma\big( \frac{|\alpha|+2+p}{p} \big) }.$ We restrict our attention to the part of $A_\alpha$ varying with respect to $\alpha$, i.e. $\binom{|\alpha|}{\alpha}\prod_{i=1}^2 \Gamma \big((\alpha_i +1)/p \big).$ Let us observe, that for $|\alpha|=4$ those coefficients attain only two different values: \begin{itemize} \item[{(1)}] $\Gamma\big(\frac{5}{p} \big)\Gamma \big(\frac{1}{p}\big)$, whenever $\alpha=(4,0)$ or $\alpha = (0,4)$. This coefficient stands by $\frac{\partial^4}{\partial x^4}$ and $\frac{\partial^4}{\partial y^4}$ in $R_4$, \item[{(2)}] $6\Gamma\big(\frac{3}{p}\big)^2$, if $\alpha=(2,2)$. This coefficient appears by $\frac{\partial^4}{\partial x^2 \partial y^2}$ in operator $R_4$. \end{itemize} Therefore $R_4$ takes a form $\Big( \frac{2}{p}\Big)^2 \Gamma \Big( \frac{|\alpha|+4}{p}\Big)^{-1} \Big([\Gamma\Big(\frac{5}{p} \Big)\Gamma \Big(\frac{1}{p}\Big) \Big( \frac{\partial^4}{\partial x^4} + \frac{\partial^4}{\partial y^4}\Big) +6\Gamma\Big(\frac{3}{p}\Big)^2 \Big( \frac{\partial^4}{\partial x^2 \partial y^2}\Big) \Big],$ which, up to a multiplicative constant, reduces to $\Delta^2= \frac{\partial^4}{\partial x^4} + \frac{\partial^4}{\partial y^4} + 2 \frac{\partial^4}{\partial x^2 \partial y^2}$ if and only if $f(p):= \frac{\Gamma\big(\frac{3}{p}\big)^2}{\Gamma\big(\frac{5}{p} \big) \Gamma \big(\frac{1}{p}\big)} = \frac{1}{3}.$ By the previous considerations this holds true for $p=2$. Let us differentiate $f$ with respect to $p$. Recall, that the formula for derivative of the gamma function stays: $\Gamma'(z) = \Gamma(z) \Big(-\frac{1}{z} -\gamma -\sum_{k=1}^\infty \frac{1}{k+z} -\frac{1}{k} \Big) = \Gamma(z) \Psi(z),$ where $\gamma$ is the Euler constant and $\Psi$ is the digamma function defined by the above equality. We use this identity to compute the following \begin{align*} f'(p) &= \frac{2\Gamma\big(\frac{3}{p}\big)^2 \Psi\big(\frac{3}{p}\big) \big(-\frac{3}{p^2}\big)\Gamma\big(\frac{5}{p}\big) \Gamma\big(\frac{1}{p}\big)}{\Gamma\big(\frac{5}{p}\big)^2 \Gamma\big(\frac{1}{p}\big)^2}\\ &\quad -\frac{\Gamma\big(\frac{3}{p}\big)^2\big[ \Gamma\big(\frac{5}{p}\big)\Psi \big(\frac{5}{p}\big)\big(-\frac{5}{p^2}\big)\Gamma\big(\frac{1}{p}\big) + \Gamma\big(\frac{5}{p}\big)\Gamma\big(\frac{1}{p}\big)\Psi\big(\frac{1}{p}\big) \big(-\frac{1}{p^2}\big) \big]}{\Gamma\big(\frac{5}{p}\big)^2\Gamma\big(\frac{1}{p}\big)^2} \\ &= f(p) \Big( -\frac{6}{p^2} \Psi\Big(\frac{3}{p}\Big) + \frac{5}{p^2} \Psi\Big(\frac{5}{p}\Big) +\frac{1}{p^2} \Psi \Big(\frac{1}{p} \Big)\Big). \end{align*} Since $f$ is positive for $p \in [1,\infty)$, we only need to investigate the sign of the second factor in the above formula: \begin{align*} &-\frac{6}{p^2} \Psi\Big(\frac{3}{p}\Big) + \frac{5}{p^2} \Psi\Big(\frac{5}{p}\Big) +\frac{1}{p^2} \Psi\Big(\frac{1}{p}\Big) \\ &= - \frac{6}{p^2}\Big( -\frac{p}{3}-\gamma - \sum_{k=1}^\infty \frac{1}{k+3/p} -\frac{1}{k}\Big)+ \frac{5}{p^2} \Big( -\frac{5}{p}-\gamma - \sum_{k=1}^\infty \frac{1}{k+5/p}-\frac{1}{k}\Big)\\ &\quad + \frac{1}{p^2} \Big( -p-\gamma -\sum_{k=1}^\infty \frac{1}{k+1/p} -\frac{1}{k}\Big)\\ &= -\frac{1}{p^2} \sum_{k=1}^\infty \frac{-6}{k+3/p} +\frac{5}{k+5/p}+\frac{1}{k+1/p}\\ &= \frac{1}{p^2}\sum_{k=1}^\infty \frac{8k}{p(k+3/p)(k+5/p)(k+1/p)} >0. \end{align*} Therefore $f$ is monotonically increasing on $[1,\infty)$ and attains value $1/3$ exactly at $p=2$. We conclude our computations with the following: $R_4 = \begin{cases} \big(\frac{2}{p}\big)^2 \Gamma \big( \frac{|\alpha|+4}{p}\big)^{-1} \Gamma\big(\frac{5}{2}\big)\Gamma\big( \frac{1}{2}\big) \Delta^2 & \text{for } p=2, \\[4pt] \big(\frac{2}{p}\big)^2 \Gamma \big( \frac{|\alpha|+4}{p}\big)^{-1} \Gamma\big(\frac{5}{p} \big)\Gamma \big(\frac{1}{p}\big) \Big( \Delta^2 + \Big( \frac{6\Gamma\big(\frac{3}{p}\big)^2} {\Gamma\left(\frac{5}{p} \right)\Gamma \big(\frac{1}{p}\big)} -2 \Big) \frac{\partial^4}{\partial x^2 \partial y^2} \Big) & \text{for } p\ne 2. \end{cases}$ We are now in a position to apply Theorem \ref{twsys}. Let $u \in \mathcal{H}(\Omega,l^p,dx)$. Then it satisfies the system of equations \eqref{system3} with $w=1$ and so \eqref{eq11} reads \begin{equation}\label{lap} \Delta u = 0, \end{equation} hence $u$ is harmonic, and its bilaplacian vanishes. Moreover $u$ has to satisfy equation of system \eqref{system3} for $j=4$, i.e.\ $R_4 u =0$. Since bilaplacian of $u$ vanishes, therefore $u$ is in fact solution to $u_{x x y y}=0$. Let us observe, that differentiating twice $\Delta u$ with respect to $x$ and $y$ respectively we obtain $u_{x x x x} + u_{x x y y} = 0 \quad \text{and} \quad u_{x x y y} + u_{y y y y} = 0.$ Therefore both $u_{x x x x}=0$ and $u_{y y y y}=0$, which means that for each fixed value of $y$ function $u(x,y)$ is a polynomial in $x$ of degree at most 3 and analogously for a fixed $x$ function $u(x,y)$ is a polynomial in $y$ of degree at most 3. Then there exist $a_i(y)$ and $b_i(x)$ for $i=0,1,2,3$ such that \begin{equation}\label{e0} u(x,y) = a_0(y) +a_1(y) x + a_2(y) x^2 +a_3(y) x^3 = b_0(x) + b_1(x) y+b_2(x) y^2 + b_3(x) y^3. \end{equation} In what follows we omit writing the arguments of $a_i$ and $b_i$. Simple calculations give us \begin{gather}\label{e1} u_{x x x x}=b_0^{(4)} +b_1^{(4)} y + b_2^{(4)} y^2 + b_3^{(4)} y^3 =0, \\ \label{e2} u_{y y y y}=a_0^{(4)} + a_1^{(4)} x +a_2^{(4)} x^2 +a_3^{(4)} x^3 =0. \end{gather} Now at each fixed $x$ in \eqref{e1} the polynomial in $y$ has to have all coefficients equal to $0$ due to the Equality of Polynomials Theorem, hence $b_i^{(4)} = 0$ for $i=1,2,3,4$. Similarly, at \eqref{e2} we set that $a_i^{(4)}=0$ for all $i=1,2,3,4$. Therefore, all of $a_i$ and $b_i$ are polynomials of degree at most $3$. Moreover, we know that $u_{x x y y}=0$. We calculate this derivative in \eqref{e0} to obtain $0=u_{x x y y} =2a_2'' + 6 x a_3'' = 2b_2'' +6yb_3 ''.$ Thus, once again we obtain that $a_i''=0$ and $b_i''=0$ for $i=2,3$, so $a_2$, $a_3$, $b_2$ and $b_3$ are in fact of degree at most $1$. By the above considerations we conclude that $u$ is a linear combination of the monomials \begin{equation}\label{e3} 1,x,y,xy,x^2,x^3,xy^2,xy^3,x^2y,x^3y,y^2,y^3, \end{equation} which solves equation \eqref{lap}. Therefore, $u$ has to be a harmonic polynomial of the form described by \eqref{e3}. The part of $u$ generated by $\{ 1,x,y,xy\}$ is already harmonic and for that reason we only need to consider $u$ being a combination of the remaining monomials in \eqref{e3}, i.e. $u=c_1 x^2 + c_2 x^3 + c_3 xy^2 +c_4 xy^3 + c_5 x^2y+c_6 x^3y + c_7 y^2 + c_8 y^3.$ Inserting $u$ to \eqref{lap} we get the following $0=2 (c_1 + c_7) + 2x(3c_2 + c_3) + 6xy(c_4+c_6) + 2y(c_5+3 c_8),$ and once again using the Equality of Polynomials Theorem we obtain that \begin{equation}\label{e33} u \in \operatorname{span} \big\{1,x,y,xy, x^2-y^2, xy^2-\frac{x^3}{3}, xy^3-x^3y,x^2y-\frac{y^3}{3} \big\}. \end{equation} Finally, let us observe that in equations of system \eqref{system3} for $j=6$ there appear only the following operators $\frac{\partial^6}{\partial x^6}, \frac{\partial^6}{\partial x^4 \partial y^2}, \frac{\partial^6}{\partial x^2 \partial y^4},\frac{\partial^6}{ \partial y^6},$ which all vanish on $u$ as in \eqref{e33}. The triviality of equations for $j>6$ follows immediately. Therefore, we summarize our discussion with the inclusion \begin{equation}\label{lp} \mathcal{H}(\Omega,l^p,dx) \subset \operatorname{span} \big\{1,x,y,xy, x^2-y^2, xy^2-\frac{x^3}{3}, xy^3-x^3y,x^2y-\frac{y^3}{3} \big\}. \end{equation} We postpone the proof of the opposite inclusion till Section \ref{51}. \subsection{The case of $l^p$ distance for $p =\infty$} To complete our illustration of Theorem \ref{twsys} we need to consider the remaining case, i.e.\ characterize functions $u$ in $\mathcal{H}(\Omega,l^p,dx)$ for $p=\infty$ by using Theorem \ref{twsys}. In this case $B(0,1)=[-1,1]^n$ in $l^\infty$ norm. Therefore, we obtain the following formula for the coefficients $A_\alpha$ in \eqref{system3}: $A_\alpha =\binom{|\alpha|}{\alpha} \int_{-1}^1 x_1^{\alpha_1} \cdot \ldots \int_{-1}^1 x_n^{\alpha_n} = \binom{|\alpha|}{\alpha}\frac{2^n}{\prod_{i=1}^n (\alpha_i +1)}.$ Then, after inserting $A_\alpha$ and dividing by the $2^n$ factor, system \eqref{system3} converts to \begin{equation} \sum_{ \substack{|\alpha|=j \\ \alpha_i \text{ even}}} \binom{|\alpha|}{\alpha} \frac{1}{(\alpha_1+1)! \cdots (\alpha_n+1) !} D^\alpha u =0. \end{equation} As in the previous subsection we restrict our attention to case $n=2$ and write out the equation for $j=2$: $\frac{1}{6} (u_{x x} + u_{y y})=0$. Hence $u$ is a harmonic function. Equation for $j=4$ is the following $\frac{1}{120} (u_{xxxx}+u_{yyyy}) + \frac{1}{6} u_{xxyy} =0,$ and can be reduced to $\Delta^2 u + 20 u_{xxyy}=0$. This, combined with an analogous discussion to the one ending the previous subsection leads us to the conclusion that inclusion \eqref{lp} holds also for $p=\infty$. \section{Converse of Theorem \ref{twsys}} Since both Theorems \ref{F-L}, \ref{twbo} and Proposition \ref{bosod} give not only the necessary, but also the sufficient condition for the mean value properties in the sense of Flatto and Bose, respectively, our next goal is to find an appropriate counterpart of these results. In case of nonconstant weights Proposition \ref{bosod} imposes an additional PDE condition on $w$, hence we expect an analogous condition. From the point of view of our further considerations, the following \emph{generalized Pizzetti formula} introduced by Zalcman in \cite{F}, will be vital. \begin{theorem}[{\cite[Theorem 1]{F}}]\label{twza} Let $\mu$ be a finite Borel measure on $\mathbb{R}^n$ with compact support and \sloppy ${\mathcal{F}(\xi)=\int_{\mathbb{R}^n}e^{-i(\xi y)}d\mu(y)}$ be the Fourier transform of the measure $\mu$. Suppose that $f$ is an analytic function on a domain $\Omega \subset \mathbb{R}^n$. Then the following equality holds \begin{equation}\label{gpiz} \int_{\mathbb{R}^n} u(x+ry) d\mu(y) = [\mathcal{F}(-rD)f](x), \end{equation} for all $x \in \Omega$ and $r>0$ such that the left-hand side exists and the right-hand side converges. The symbol~$D$ is given by $D:=-i \big( \frac{\partial}{\partial x_1},\ldots,\frac{\partial}{\partial x_n} \big)$. \end{theorem} \begin{remark} \rm Formula \eqref{gpiz} is the main tool used in the proof of Theorem \ref{twod}, hence we need to assume analyticity of weight $w$. From a result by \L ysik \cite{GL1} the Pizzetti formula on Euclidean balls is valid exactly for analytic functions. Therefore, dropping the analyticity assumption of $w$ would require finding a different proof of Theorem \ref{twod}. \end{remark} Theorem \ref{F-L} is a special case of Theorems \ref{twsys} and \ref{twod} for $w\equiv 1$. Proposition~\ref{bosod} is generalized by Theorem \ref{twod} due to the following lemma. \begin{lemma}\label{lembo} Suppose that $w \in C^{2l}(\Omega)$ solves the equation \begin{equation}\label{gef} \Delta^l w + a_{l-1} \Delta^{l-1} w + \ldots + a_1 \Delta w + a_0 w = 0, \end{equation} where all $a_i \in \mathbb{C}$ for $i=0,1,\ldots,l-1$. Then $w$ is analytic in $\Omega$. \end{lemma} \begin{proof} We prove the lemma by the mathematical induction with respect to $l$. Recall the following fact (see \cite[p. 57]{fryc}): Suppose that $w \in C^2(\Omega)$ solves the equation \begin{equation}\label{daur} L w + \lambda w = \varphi, \end{equation} where $L$ is elliptic with analytic coefficients and $\varphi$ is analytic in $\Omega$, $\lambda \in \mathbb{C}$. Then $w$ is analytic in $\Omega$. If $l=1$, then we use the above regularity fact with $L=\Delta$, $a_0=\lambda$ and $\varphi \equiv 0$. Now let us assume that the assertion holds for $l-1$ and consider $w$ as in \eqref{gef}. By adding and subtracting the appropriate terms we may rewrite this equation as follows with any $\lambda \in \mathbb{C}$ and given $a_0,a_1,\ldots,a_{l-1} \in \mathbb{C}$: \begin{align*} 0&=\Delta^{l-1} (\Delta w + \lambda w) +(a_{l-1}-\lambda) \Delta^{l-2} (\Delta w + \lambda w ) \\ &\quad + (a_{l-2}-\lambda(a_{l-1}-\lambda)) \Delta^{l-3}(\Delta w + \lambda w) + \ldots\\ &\quad + (a_1 - \lambda(a_2-\lambda(\ldots) ) ) (\Delta w + \lambda w) + \big(a_0 -\lambda(a_1 - \lambda(a_2-\lambda(\ldots) ) )\big) w. \end{align*} Since the factor in the last $w$-term is a complex polynomial in $\lambda$, one can choose such $\lambda$, so that this last factor standing by $w$ in the equation vanishes (e.g. take $\lambda$ to be one of the roots). We use the assumption for $l-1$ to obtain that $\Delta w + \lambda w$ is an analytic function, denoted by $\varphi$, i.e.\ $\Delta w + \lambda w = \varphi$. This observation together with the regularity observation allow us to conclude the proof. \end{proof} Before proving Theorem \ref{twod} we need to show the following regularity result for strongly harmonic functions. \begin{lemma}\label{lem4} Suppose that $\Omega \subset \mathbb{R}^n$ is open, $w$ is a positive analytic function and $d$ is induced by a norm. Then any $u \in \mathcal{H}(\Omega,d,wdx)$ is analytic as well. \end{lemma} \begin{proof} By Theorem \ref{twsys} function $u$ is a weak solution to the equation for $j=2$ of system \eqref{system3}. Let us show that this equation is strongly elliptic. Take $\xi \in \mathbb{R}^n$ and consider the second order terms of this equation. By \eqref{eq14} we obtain for all $y \in \Omega$ that \begin{align*} \sum_{|\alpha|=2} A_\alpha w(y) \xi^\alpha =w(y) \int_{B(0,1)} (x \cdot \xi)^2 dx \geq w(y) \int_{\|x\|_2 \leq \varepsilon } (x \cdot \xi)^2 dx, \end{align*} where the last estimate holds with some $\varepsilon>0$ since $d$ is equivalent to the Euclidean distance. Next, observe that \begin{align*} \int_{\|x\|_2 \leq \varepsilon } (x \cdot \xi)^2 dx &= \int_{\|x\|_2 \leq \varepsilon} \cos^2 \angle( x,\xi) \|x\|_2^2 \|\xi\|_2^2 dx \\ &= \|\xi\|_2^2 \int_{\|x\|_2 \leq \varepsilon} \cos^2 \angle( x,\xi) \|x\|_2^2 dx \\ &= \theta \|\xi\|_2^2, \end{align*} where $\theta>0$ is defined by the above equality. It does not depend on $\xi$ because of the symmetry of the Euclidean ball. Hence $\sum_{|\alpha|=2} A_\alpha w(y) \xi^\alpha \geq \theta w(y) \|\xi\|_2^2, \quad y \in \Omega.$ Therefore the operator $\sum_{|\alpha|=2} A_\alpha w D^\alpha$ is strongly elliptic. We apply the regularity result \eqref{daur} to obtain that $u$ is analytic. \end{proof} \begin{proof}[Proof of Theorem \ref{twod}] We need to show the equality \begin{equation}\label{equ8} u(x) \int_{B(0,1)} w(x+ry) dy = \int_{B(0,1)} u(x+ry)w(x+ry) dy. \end{equation} To prove it we use the generalized Pizzetti formula for a measure $\mu$ being the normalized Lebesgue measure on the unit ball. Then \begin{align*} \mathcal{F}(\xi)& = \int_{B(0,1)} e^{-i \xi y} dy = \sum_{j=0}^\infty \frac{(-i)^j}{j!} \int_{B(0,1)} (\xi y )^j dy \\ &= \sum_{j=0}^\infty \frac{(-i)^j}{j!} \sum_{|\alpha|= j }A_\alpha \xi^\alpha = \sum_{\alpha \in \mathbb{N}^n} \frac{(-i)^{|\alpha|}}{|\alpha|!} A_\alpha \xi^\alpha, \end{align*} where $A_\alpha =\binom{|\alpha|}{\alpha} \int_{B(0,1)} y^\alpha dy$. We apply Theorem \ref{twza} twice: to $w$ and $uw$ to obtain \begin{gather} \label{e43} \int_{B(0,1)} w(x+ry) dy = \sum_{\alpha \in \mathbb{N}^n} \frac{r^{|\alpha|}}{|\alpha|!} A_\alpha D^\alpha w(x),\\ \label{e44} \int_{B(0,1)} u(x+ry) w(x+ry) dy = \sum_{\alpha \in \mathbb{N}^n} \frac{r^{|\alpha|}}{|\alpha|!} A_\alpha D^\alpha (u(x)w(x)), \end{gather} Multiply \eqref{e43} by $u(x)$ and subtract from it \eqref{e44} to obtain \begin{align*} & u(x) \int_{B(0,1)} w(x+ry) dy-\int_{B(0,1)} u(x+ry) w(x+ry) dy \\ &= \sum_{\alpha \in \mathbb{N}^n} \frac{r^{|\alpha|}}{|\alpha|!} A_\alpha \big( u(x) D^\alpha(w(x)) - D^\alpha (u(x)w(x)) \big) \\ &= \sum_{j=0}^\infty \frac{r^j}{j!} \sum_{|\alpha|=j} A_\alpha \big( u(x) D^\alpha(w(x)) - D^\alpha (u(x)w(x)) \big) =0, \end{align*} where in the last step we appeal to \eqref{system3}. Thus, the proof is complete. \end{proof} \subsection{Applications of Theorem \ref{twod}}\label{51} In the last part of our work we present some of the consequences of Theorem \ref{twod}. Consider the set $\Omega \subset \mathbb{R}^2$, metric $d$ induced by the $l^p$ norm for $1\leq p \leq \infty$, $p \ne 2$ and $\mu$ being the Lebesgue measure. From the computations summarized in \eqref{e3} we know, that the set of solutions to system~\eqref{system3} is $\operatorname{span} \big\{1,x,y,xy, x^2-y^2, xy^2-\frac{x^3}{3}, xy^3-x^3y,x^2y -\frac{y^3}{3} \big\}$. This together with Theorem \ref{twod} allows us to augment observation \eqref{lp} with the converse inclusion: $\mathcal{H}(\Omega, l^p,dx) =\operatorname{span} \big\{1,x,y,xy, x^2-y^2, xy^2-\frac{x^3}{3}, xy^3-x^3y,x^2y-\frac{y^3}{3} \big\}.$ Notice, that the dimension of $\mathcal{H}(\Omega, l^p,dx)$ is equal to 8. As mentioned in the introduction, in $\mathbb{R}^3$ \L ysik \cite{GL} computed $\dim \mathcal{H}(\Omega, l^p,dx)= 48$. Those numbers coincide with $2^n n!$ - the number of linear isometries of $(\mathbb{R}^n,l^p)$, which is discovered in \cite{LI}. We believe that there is a link between the dimension of the space $\mathcal{H}(\Omega, d,\mu)$ and the number of linear isometries, still to be examined. Moreover, let us consider the metric measure space $(\Omega,l^2,wdx)$ for the analytic weight function $w$. Then, by Theorems \ref{twsys} and \ref{twod} and Section \ref{ldwa} we know that $\mathcal{H}(\Omega,l^2,wdx )$ consists exactly of solutions to the system of equations \begin{equation}\label{e5} \Delta u \Delta^j w +2 \nabla u \nabla (\Delta^j w)=0, \quad \text{for } j=0,1,\ldots. \end{equation} Let us observe, that $u$ solves also infinitely many other systems of equations, obtained from \eqref{e5} by excluding $l \in \mathbb{N}$ initial equations $\Delta u \Delta^{j+l} w +2 \nabla u \nabla (\Delta^{j+l} w) = \Delta u \Delta^j(\Delta^l w) +2 \nabla u \nabla (\Delta^j(\Delta^l w))=0,$ for $j=0,1,\ldots$. Therefpre, $u$ is strongly harmonic in countably many metric measure spaces $(\Omega,l^2,\Delta^l w dx)$ for all $l \in \mathbb{N}$. In other words, function $u$ has infinitely many mean value properties, with respect to different weighted Lebesgue measures $d\mu = \Delta^l w dx$ for all $l \in \mathbb{N}$, whenever $\Delta^l w$ are positive. \subsection*{Acknowledgements} The author is particularly grateful to his academic advisor Tomasz Adamowicz for important comments, fruitful discussions and valuable lessons in many cases extending results of this article. \begin{thebibliography}{00} \bibitem{H} T. Adamowicz, M. Gaczkowski, P. G\'orka; \textit{Harmonic functions on metric measure spaces.} Rev. Math. Complut., 32, (2019), 141--186. \bibitem{H1} T. Adamowicz, B. Warhurst; \textit{Mean value property and harmonicity on Carnot-Carathéodory groups.}, Potential Anal., (2018), doi: 10.1007/s11118-018-9740-4. \bibitem{LL2} A. Arroyo, J. G. Llorente; \textit{A priori Hölder and Lipschitz regularity for generalized p-harmonious functions in metric measure spaces.} Nonlinear Anal., 168, (2018), 32--49. \bibitem{LL3} A. Arroyo, J. G. Llorente; \textit{On the Dirichlet problem for solutions of a restricted nonlinear mean value property.} Differential Integral Equations, 29, (2016), no. 1-2, 151--166. \bibitem{B2} A. K. Bose; \textit{Functions satisfying a weighted average property.} Trans. Amer. Math. Soc., 118, (1965), 472--487. \bibitem{B3} A. K. Bose; \textit{Functions satisfying a weighted average property II.} Trans. Amer. Math. Soc., 124, (1966), 540--551. \bibitem{B4} A. K. Bose; \textit{Generalized eigenfunctions of the Laplace operator and weighted average property.} Proc. Amer. Math. Soc., 19, (1968), 55--62. \bibitem{LI} D. Andrica, V. Bulgarean; \textit{Some remarks on the group of isometries of a metric space.} Nonlinear analysis, 57--64, Springer Optim. Appl., 68, Springer, New York, (2012). \bibitem{FAA} P. Duarte, M. J. Torres; \textit{A discrete Faà di Bruno's formula.} Contrib. Discrete Math., 7, (2012), no. 2, 72--83. \bibitem{D1} J. Edwards; \textit{A Treatise on the Integral Calculus. Vol. II} Chelsea Publishing Company, New York, (1922). \bibitem{E2} L. C. Evans; \textit{Partial differential equations.} Graduate Studies in Mathematics 19, American Mathematical Society, Providence, RI, (1998). \bibitem{EG} L. C. Evans, R. F. Gariepy; \textit{Measure theory and fine properties of functions.} Studies in Advanced Mathematics., CRC Press, Boca Raton, FL, (1992). \bibitem{FLM} F. Ferrari, Q. Liu, J. J. Manfredi; \textit{On the characterization of $p$-harmonic functions on the Heisenberg group by mean value properties.} Discrete Contin. Dyn. Syst., 34, (2014), no. 7, 2779--2793. \bibitem{ferrari} F. Ferrari, A. Pinamonti; \textit{Characterization by asymptotic mean formulas of q-harmonic functions in Carnot groups.} Potential Anal., 42, (2015), no. 1, 203--227. \bibitem{F2} L. Flatto; \textit{Functions with a mean value property.} J. Math. Mech., 10, (1961), 11--18. \bibitem{F3} L. Flatto; \textit{Functions with a mean value property II.} Amer. J. Math., 85, (1963), 248--270. \bibitem{F1} L. Flatto; \textit{The converse of Gauss's theorem for harmonic functions.} J. Differential Equations, 1, (1965), 483--490. \bibitem{F5} A. Friedman, W. Littman; \textit{Functions satisfying the mean value property.} Trans. Amer. Math. Soc., 102, (1962), 167--180. \bibitem{J} M. Gaczkowski, P. G\'orka; \textit{Harmonic Functions on Metric Measure Spaces: Convergence and Compactness.} Potential Anal., 31, 203--214, (2009). \bibitem{A} W. Hansen, N. Nadirashvili; \textit{A converse to the mean value theorem for harmonic functions.} Acta Math., 171, (1993), no. 2, 139--163. \bibitem{fryc} F. John; \textit{Plane waves and spherical means applied to partial differential equations.} Interscience Publishers, New York-London, (1955). \bibitem{LL1} J. G. Llorente; \textit{Mean value properties and unique continuation.} Commun. Pure Appl. Anal., 14, (2015), no. 1, 185--199. \bibitem{GL1} G. \L ysik; \textit{A characterization of real analytic functions.} Ann. Acad. Sci. Fenn. Math., 43, (2018), no. 1, 475--482. \bibitem{GL} G. \L ysik; \textit{personal communication.} (2018). \bibitem{man} J. J. Manfredi, M. Parvianen, J. D. Rossi; \textit{On the definition and properties of p-harmonious functions.} Ann. Sc. Norm. Super. Pisa Cl. Sc., 11, (2013), 215--241. \bibitem{OO} R. Schneider; \textit{Convex bodies: the Brunn-Minkowski theory. Second expanded edition.} Encyclopedia of Mathematics and its Applications, 151. Cambridge University Press, Cambridge, (2014). \bibitem{DS} D. Schymura; \textit{An upper bound on the volume of the symmetric difference of a body and a congruent copy.} Adv. Geom., 14, (2014), no. 2, 287--298. \bibitem{F} L. Zalcman; \textit{Mean values and differential equations.} Israel J. Math., 14, (1973), 339--352. \end{thebibliography} \end{document}