Abstract: We introduce a new boundary Harnack principle in Lipschitz domains for equations with a right hand side. Our approach, which uses comparisons and blow-ups, seems likely to adapt to more general domains, as well as other types of operators.

We consider operators of the form $$ \mathcal{L} u = (a^{ij}u_i)_j + b^i u_i + cu , $$ with the following ellipticity conditions $$ \Lambda^{-1} |\xi|^2 \leq \langle a^{ij}(x) \xi, \xi \rangle \leq \Lambda |\xi|^2, $$ for some $\Lambda >0$ and for all nonzero $\xi \in \R^n$. Furthermore, $a^{ij}(x)$ is a real $n \times n$ matrix. For the lower order terms we assume $$|c(x)|, \sum |b^i(x)| \leq \Lambda-1$$ and that $c(x) \leq 0$. We say $u \in \mathcal{S_{\mathcal{L}}}(\Omega_{L,R})$ if $$ \begin{aligned} \mathcal{L} u(x) &=0 \text{ in } \Omega_{L,R} \cap B_R ,\\ u(x)&=0 \text{ on } \Omega_{L,R}^c\cap B_R, \end{aligned} $$ and that $u \in \mathcal{S_{\mathcal{L}}}(\Omega_{L,R},d^{\gamma})$ if $$ \begin{aligned} |\mathcal{L} u(x)| &\leq (\text{dist}(x,\partial \Omega_{L,R} \cap B_R))^{\gamma} \ \text{ in } \Omega_{L,R} , \\ u(x)&=0 \ \text{ on } \Omega_{L,R}^c\cap B_R. \end{aligned} $$ To apply the Hölder continuity estimates for elliptic operators we will require that $\gamma>-2/n $. Since the boundary is Lipschitz, this will ensure the correct integrability assumptions for the right hand side.

Theorem:

Let $0 \in \partial \Omega_L$ with $L< M$, and fix $x^0 \in \Omega_L$. Assume further that $B_1 \cap \{x_n > 1/4\} \subseteq \Omega_L$. Assume $u,v \geq 0$ and $u,v \in \mathcal{S}(\Omega_{L},d^{\gamma})$ with $\mathcal L u, \mathcal L v \leq 0$ and $u(x^0)=v(x^0)=1$, and also assume that $2-\alpha+\gamma>0$, with $\gamma>-2/n$ and $\alpha$ be such that $$ \sup_{B_{r}(x)} u \geq c_1 u(e_n/2) r^{\alpha}. $$

Then there exists a uniform constant $C>0$ (depending only on dimension $n$, Lipschitz constant $M$, $M-L$, and dist$(x^0,\partial \Omega_L)$) such that \begin{equation} \label{e:ine3} C^{-1} v(x) \leq u(x) \leq Cv(x) \end{equation} for all $x \in B_{1/2}$.

The aim of this paper is to establish a link between quadrature surfaces (potential theoretic concept) and sandpile dynamics (Laplacian growth models). Given a measure $ \mu $ in $\mathbb{R}^d$ and a domain $\Omega$ containing the support of $\mu$, we call $\partial \Omega$ a quadrature surface for $\mu $ if one has $ \int_{\partial \Omega } h(x) d\mathcal{H}^{d-1} (x) = \int h(x) d\mu (x) $ for all functions $ h $ harmonic on $ \overline{\Omega} $ [Sh1] and [Sh2] by H. Shahgholian). For example, if $\mu$ is a Dirac mass concentrated at some point $x\in \mathbb{R}^d$, then a sphere with center at $x$ and of an appropriate radius would give a quadrature surface. There is no complete description of measures $\mu$ admitting a quadrature surface identity, but in a special case of a finite sum of Dirac masses, for instance, the existence is due to Shahgholian.

In this paper we show that one can obtain a quadrature surface as a scaling limit of a certain growth model on the lattice $\mathbb{Z}^d $, where $ d\geq 2 $. Our motivation comes from an already established connection between quadrature domains (these are domains admitting a similar integral identity as above, but with integration against the Lebesgue measure over the domain instead of the surface) and Laplacian growth models by Levine and Peres (see [LP1] and [LP2] ). They introduce a model called divisible sandpile, and prove that for initial distribution of finitely many points its scaling limit is a quadrature surface ( a smash sum of Euclidean balls, precisely, a concept after Diaconis and Fulton) A shift from a domain to a lower-dimensional object, a hypersurface in our case, seems to necessitate an entirely different approach.

To define our model, start with a distribution of mass $\mu_0 $ on $ \mathbb{Z}^d $, i.e. a bounded non-negative function of finite support, let $n = \sum_{x\in \mathbb{Z}^d} \mu_0(x) $ be its total mass, and fix a threshold $ m>0 $. One by one we pick a vertex $x \in \mathbb{Z}^d$ and if it carries mass larger than $m$, distribute the excess of mass from $m$ evenly among the $2d$ lattice neighbors of $x$. Otherwise, if the mass carried by $x$ is bounded by $m$ but the cumulative emissions of mass from $x$ prior to a given epoch of time exceed $\frac 1m n^{2/d}$ distribute the entire mass of $x$ evenly among its lattice neighbors. If none of the two conditions hold for $x$ we move on to another vertex. This mass redistribution procedure is called toppling of $x$. We show that provided any vertex of $\mathbb{Z}^d$ is chosen infinitely many times, the mass redistribution process will reach a stable state after countably many topplings, and regardless of the order of the topplings the final configuration will be the same (Abelian property). Informally, one may think about the proposed model, as a deformation of a divisible sandpile model of Levine and Peres, where we force the mass to accumulate on a $\frac 1m $-sub-level sets of the odometer. This model is partially inspired by singular perturbation problems in PDEs which are known to converge to Bernoulli type free boundaries.

The key concept in our analysis is the odometer function $u:\mathbb{Z}^d \rightarrow \mathbb{R}_+$ where for each $x\in \mathbb{Z}^d$ the value $u(x)$ shows total emissions of mass from $x$ during the lifetime of the process. In terms of this function, the discrete PDE problem solved by the model reads $$\Delta u(x) + \mu_0(x) \leq m \ \ \ \text{ everywhere on } \mathbb{Z}^d $$ and $$\Delta u(x) = m \mathbb{I}_{ \{ 0 < u < \frac{1}{m} n^{2/d} \} } (x) - \mu_0(x) \ \ \ \text{ for all } x \in \{u > 0\},$$ where $\Delta$ is the discrete (normalized) Laplacian. We show that the odometer $u$ is precisely the smallest super-solution to the above problem. This characterization combined with combinatorial moving plane techniques, which we introduced in our earlier paper, allows to get a certain discrete monotonicity properties for the odometer function for a single source initial distributions, using which we show that the set of visited sites of the model in $\mathbb{Z}^d$ grows proportional to $n^{1/d}$, and the entire mass of the system is being distributed to an annular ring of thickness $\sim 1/m$.

We prove the existence of the scaling limit of the model generated by a single source, and fixed threshold $m$. As we show, this limit is a ball, with the entire mass of the system being redistributed onto a narrow annular ring of thickness $\frac 1m$ near the boundary of the ball. With compactness arguments, we conclude that there is also a scaling limit for a subsequence of the odometers, if the threshold $m$ tends to infinity along with the scale of the model. That limit is spherical, with the entire mass of the system being uniformly redistributed onto the boundary of that ball. This gives a quadrature surface for a single source mass distributions.

We introduce a new lattice growth model, which we name boundary sandpile and which amounts to potential theoretic redistribution of the given continuous mass on the lattice $\mathbb{Z}^d$ onto the combinatorial (free) boundary of some unknown domain, which is being determined by the dynamics of the model.

Our initial motivation was triggered by an intriguing connection established by Levine and Peres, between quadrature domains, a potential theoretic concept, and divisible sandpile, a growth model on the lattice. In this paper we attempted to create a sandpile dynamics (a growth model on $\mathbb{Z}^d$) which would have a quadrature surface as its scaling limit. On the way to this quest, which is yet to be fulfilled, we introduce a new growth model, the boundary sandpile, which seems to represent a new alluring phenomenon not encountered previously.

To define boundary sandpile model, assume we have concentration of (continuous) mass $n > 0$ at the origin of $\mathbb{Z}^d$ $(d\geq 2)$. Then, at any discrete time $k\geq 0$, we say that a site $x\in \mathbb{Z}^d$ is visited, if at a certain time before (and including) $k$ it carried positive amount of mass. Next, at time $k$ we call a site $x$ unstable, if either $x$ is in the interior of the visited sites and has positive mass, or $x$ is on the (discrete) boundary of the set of visited sites, and has mass larger than $n^{1/d}$ (this threshold is due to scaling considerations). Then, any unstable site can topple, by evenly distributing all its mass among its neighboring vertices on the lattice. The process terminates, if there are no unstable sites. Clearly, this redistribution never stops in finite time, except for trivial cases, and it is not clear a priori, if the order in which one topples the unstable sites, can affect in any way the limiting growth cluster. What we showed in the paper, is that in the limit, i.e. as the discrete time goes to infinity, the growth process finds a stable configuration provided all unstable sites are toppled infinitely often; moreover, the limiting configuration is independent of the toppling sequence, and in this sense the model is abelian. By definition, the total mass of the system is being redistributed to a combinatorial free boundary determined by the underlying dynamical system. Some limiting shapes of the sandpile are illustrated below.

Absorbing part of the topplings into linear operators acting on graphs, we show that the model reaches a stable state in finitely many steps. We also prove that the final set of visited sites is the smallest domain (with respect to inclusions) in $\mathbb{Z}^d$ where the process can stabilize itself. The analysis is carried out in terms of a special function $u$, called the odometer, which is defined on $\mathbb{Z}^d$ and for each $x\in \mathbb{Z}^d$, $u(x)$ shows the total mass emitted from $x$ during the entire life-time of the process. We show a uniform (with respect to initial mass $n$) Lipschitz estimate on the odometer function, which (in part) implies Lipschitz regularity of the free boundary of the sandpile.

The underlying ideas of the analysis are largely motivated by potential theory, free boundary problems, and elliptic PDEs, building on top of these their combinatorial and probabilistic (random walk) counterparts. Our methods of this paper had already shown to be useful in other situations too. For instance, combined with the Least Action Principle for sandpiles proved by Fey, Levine, and Peres, the techniques developed for boundary sandpile enabled us to show for the first time, that the boundary of the scaling limit of the classical Abelian sandpile model with a single source, is Lipschitz regular in dimensions $d\geq 2$. To prove this, we show that the Least Action Principle combined with the discrete variant (considered in our paper) of the celebrated moving plane method of A.D. Aleksandrov, implies monotonicity of the discrete odometer function of the sandpile with respect to special directions determined by mirror symmetries of the unit cube of $\mathbb{R}^d$. We then observe, that this monotonicity of the odometer is being transferred to the scaling limit of the model, the existence of which was established only recently by Pegden and Smart. Relying on directional monotonicity of the odometer function, we show that at each point of the free boundary of the scaling limit, there are $d$-linearly independent directions along which the odometer of the Abelian sandpile is monotone. The latter allows us to show that at each point of the free boundary there is a double-cone (of uniform size) with the vertex at that point such that the free boundary remains entirely outside the interior of the cone. That condition implies Lipschitz regularity.

We close the paper with some open problems.

We study the local behavior of the nodal sets of the solutions to elliptic quasilinear equations with nonlinear conductivity part, $\hbox{div}(A_s (x, u)u) = \hbox{div} { \mathbf f (x)}$, where $A_s (x, u)$ has " broken " derivatives of order $s \geq 0$, such as $A_s (x, u) = a(x) + b(x)(u +)^s$ , with $(u +)^0$ being understood as the characteristic function on $ \{u > 0\}$. The vector $\mathbf f (x)$ is assumed to be $C^\alpha$ in case $s = 0$, and $C^{ 1,\alpha}$ (or higher) in case $s > 0$. Using geometric methods, we prove almost complete results (in analogy with standard PDEs) concerning the behavior of the nodal sets. More exactly, we show that the nodal sets, where solutions have (linear) nondegeneracy, are locally smooth graphs. Degenerate points are shown to have structures that follow the lines of arguments as that of the nodal sets for harmonic functions, and general PDEs.

We study the minimum problem for the functional \begin{equation*} \int_{\Omega}\bigl( \vert \nabla {\mathbf{u}} \vert^{2} + Q^{2}\chi_{\{\vert {\mathbf{u}} \vert>0\}} \bigr)dx \end{equation*} with the constraint $u_i\geq 0$ for $i=1,\cdots,m$ where $\Omega\subset\mathbb{R}^{n}$ is a bounded domain and ${\mathbf{u}}=(u_1,\cdots,u_m)\in H^{1}(\Omega;\mathbb{R}^{m})$.

Using an array of technical tools, from geometric analysis for the free boundaries, we reduce the problem to its scalar counterpart and hence conclude similar results as that of scalar problem. This can also be seen as the most novel part of the paper, that possibly can lead to further developments of free boundary regularity for systems.