\section{}

Starting from (6), we wish to calculate the minimum action for a system governed by the effective lagrangian:

\begin{equation}
L_{FW}=\frac{1}{2}\left(\frac{\dot{x}-f(x)}{g(x)}\right)^2.
\end{equation}

By doing a legendre transformation we can go to the associated hamiltonian form:

\begin{equation}
H_{FW}(p,x)=p\dot{x}-L_{FW}=p\left(p\frac{g^2(x)}{2}+f(x)\right).
\end{equation}

Since the lagrangian does not depend explicitely on time, the Hamilton-Jacobi equation for the problem reduces to the simple form $H_{FW}(\partial_xS_{FW},x)=0$, where $S_{FW}$ is the action for the problem. This leads to the differential equation for $S_{FW}$:

\begin{equation}
\partial_xS_{FW}=-2\frac{f(x)}{g^2(x)},
\end{equation}
shown in equation (7) of the text.
\end{document}


\section{Introduction}
Among the numerous discoveries of the \textit{Kepler}{} mission, the detection of circumbinary (CB) planets Kepler 16 through the multiple system Kepler 47 \citep{Doyle11, Welsh12, Orosz12,  Orosz12a, Schwamb12} has  propelled a renewed theoretical effort in explaining how such planets would be assembled in a binary environment. In presence of a close binary companion, protoplanetary disks can become a rather hostile planetary nursery in both circumstellar and circumbinary configurations. The difficulties encountered by the standard core accretion paradigm in the binary environment are several, including truncation, mass loss and relatively fast dispersal of disks in close binaries \citep[e.g.][]{Duchene10, Kraus12}, possible vaporization of grains in dynamically excited disks \citep{Nelson00} and impaired planetesimal growth into protoplanets \citep{Marzari00, Moriwaki04, Thebault04, Thebault06, Scholl07, Paardekooper08, Thebault11}.

The latter ``planetesimal bottleneck'' is a robust consequence of the interplay between the gravitational perturbations of the stellar companion (which stirs the planetesimal disk and acts to raise eccentricities) and the aerodynamic drag from a putative protoplanetary disk. Since aerodynamic drag tends to both damp planetesimal eccentricities and align planetesimal orbits in a size-dependent fashion, planetesimals of different sizes will collide on different phases, resulting in high collisional speeds which lead to destructive (rather than accreting) events.
Recently, \citet{Meschiari12} (hereafter M12) investigated planetesimal accretion in the Kepler-16 system, using $N$-body simulations (coupled with drag from a static background disk) which track planetesimal collisions throughout a range of semi-major axes over $10^5$ years. The census of planetesimal collisions indicated that regions inside 4 AU ($\approx 20 a_B$, where $a_B$ is the semi-major axis of the central binary) were dominated by destructive events, and therefore hostile to planet formation. Consequently, we posited that Kepler 16 could have plausibly assembled outside the forbidden region and subsequently migrated inwards through tidal interaction with the protoplanetary disk, later stopping close to the inner edge of the disk. This scenario is supported by the evolution of planetary cores in hydrodynamical simulations \citep{Pierens07, Pierens08}, which also suggested Jupiter-mass planets in CB configurations should be rare (in accordance with the observed \textit{Kepler}{} sample). \citet{Paardekooper12} investigated planet formation in the Kepler 16, 34 and 35 systems using a similar approach, additionally including self-consistent planetesimal formation and destruction. They reached analogous conclusions, and asserted that \emph{in-situ} formation was unlikely, even under the most favorable conditions.

The approach of M12, similarly to previous investigations, neglected several physical ingredients for the sake of computational expediency. Indeed, fully self-consistent simulations which include the hydrodynamical response of the protoplanetary disk have shown that additional oscillations in eccentricity and longitude of pericenter of the planetesimals might be introduced from the development of bulk eccentricity and spiral perturbations in the disk. However, their magnitude might depend somewhat on the details of the hydrodynamical simulation \citep[e.g.][]{Paardekooper08, Marzari12, Muller12}.

In this Letter, we consider magnetohydrodynamical turbulence driven by the magnetorotational instability (MRI) as an additional source of perturbations on the planetesimal disk. MRI-driven turbulence \citep{Balbus91} is thought to be the likely source of anomalous viscosity in protoplanetary disks \citep[e.g.][]{Armitage98}, and influence how planetesimal formation \citep[e.g.][]{Johansen07}, planetesimal accretion \citep[e.g.][]{Ogihara07, Ida08, Nelson10}, and planetary migration \citep[e.g.][]{Nelson03, Laughlin04, Nelson04, Baruteau10} proceed. Previous studies simulating the dynamics of planetesimals embedded in turbulent disks in single-star environments showed that in fully MRI-active disks, the velocity dispersion of the planetesimals is significantly raised by gravitational perturbations induced by density fluctuations. For typical disk parameters and nominal turbulence strength, km-sized planetesimals might be in a highly erosive regime. Such vigorous turbulence, however, might not be appropriate to the midplane of realistic protoplanetary disks, which is thought to be dominated by a ``dead zone'' with near-laminar flow \citep{Gammie96}. \citet{Gressel11} (hereafter G11) presented the results of  stratified, 3D MHD simulations (including a substantial dead zone); planetesimals embedded in the disk midplane at  5 AU experienced a significantly reduced excitation of their eccentricities (by a factor $\approx$ 10-20). Therefore, they concluded dead zone represent ``safe havens'' for the growth of km-sized planetesimals.

Although the reduced amplitude of the random velocities excited by turbulent fluctuations is potentially conducive to planetesimal accretion in the single-star environment, the situation is more complicated in the binary environment.  In the latter case, a very precise alignment between the planetesimal orbits is crucial to attaining low encounter speeds despite the substantial eccentricity of the planetesimal orbits. Random kicks diffuse planetesimals out of alignment; although aerodynamic drag would attempt to restore alignment, it will do so on a timescale that is size-dependent, once again differentially phasing planetesimals of different sizes. Finally, high-frequency radial oscillations of the planetesimal eccentricity (as would be caused by stochastic torques) would potentially lead to further orbital crossing, leading to high encounter velocities. Therefore, we anticipate that planet formation might be strongly perturbed, or even inhibited, despite the smaller turbulent amplitudes appropriate to the midplane dead zone. We couple the $N$-body code of M12 with a numerical model that approximates the stochastic torquing arising from MRI turbulence. Our analysis shows that this additional source of perturbations is potentially damaging to planet formation in the outer disk, which was previously thought to be relatively protected from destructive impacts and therefore a plausible location for core assembly \citep[M12;][]{Paardekooper12}.

The plan of this Letter is as follows. In \S \ref{sec:setup}, we briefly discuss our numerical model for the gravitational torques arising from MHD turbulence. In \S \ref{sec:results} we show the results of our simulations, and discuss them in the context of planet formation in \S \ref{sec:conc}.

\section{Numerical setup}\label{sec:setup}
In this Letter, we study planet formation by sampling planetesimal collision events between 4 and 10 AU, and compare them to velocity thresholds corresponding to destructive impacts. We consider the orbital elements of the Kepler-16 system as our prototypical configuration, and to facilitate comparison with previous investigations \citep[M12;][]{Paardekooper12}. We refer the reader to M12 for a description of the numerical code, initial setup and velocity thresholds.

\subsection{Turbulent model}
\begin{figure}
\epsscale{1}
\plotone{f1.eps}
\caption{Sample growth plot of the eccentricity dispersion for models A, B1-4 (grey lines, from top to bottom). The dashed lines indicate the best-fit eccentricity dispersion measured in the full MHD simulations of G11, models A, D1 and D2 (from top to bottom).}\label{fig:turb_ecc}
\end{figure}
We use the analytical prescription of \citet{Laughlin04} (hereafter L04) to model torques arising from density fluctuations in the disk; we include corrections to the formulation from \citet{Ogihara07} \citep[also used in ][]{Ida08, Baruteau10}, and introduce small modifications that reduce its computational cost. In this model, density fluctuations are forced by a potential $\Phi$, a sum of turbulent $m$-fold modes
\begin{equation}
\Phi = \gamma r^2 \Omega^2 \sum_i \xi_i R(r, r_{c, i}, \sigma_i)\ T(\tilde{t}_i)\ \cos(m\varphi - \varphi_{c, i} - \Omega_{c, i} \tilde{t}_i)\label{eqn:turb}
\end{equation}
where $\gamma$ sets the overall turbulent amplitude, $r_{c, i}$, $\varphi_{c, i}$ and $\Omega_{c, i}$ are the radial center, phase and angular velocity of the mode (picked randomly in the disk), $m$ is sampled from a lograndom distribution between 1 and 6, $\sigma_i = \pi r_{c, i} / 4m$ is the radial extent of the mode, $\xi_i$ is sampled from a Gaussian distribution of unit variance, $\tilde{t}_i$ is the lifetime of the mode normalized by a timescale $\Delta t_i$ and $R$ and $T$ are two Gaussian-like functions centered around $r_{c, i}$ and $0.5 \Delta t_i$, respectively. Each mode has a limited lifetime $\Delta t_i$; following \citet{Baruteau10}, we reduce the lifetimes by a factor of 10 from the prescription of L04 in order to better match the autocorrelation timescale of 3D MHD simulations.

The actual force on the planetesimals will arise from the gravitational force of the density fluctuations induced by Equation \refp{eqn:turb}.  In order to proceed without the full hydrodynamical machinery, L04 used a WKB analysis to derive the following scaling for the RMS torque on a planetesimal of mass $\mathcal{M}_{pl}$:
\begin{equation}
\tau_T = C \gamma r^4 \Omega^2\Sigma_g \frac{\mathcal{M}_{pl}}{\mathcal{M}_*}\ \label{eqn:torque}
\end{equation}
where $\mathcal{M}_*$ is the mass of the central object (in our case, the total binary mass) and $\Sigma_g$ is the unperturbed surface density of the protoplanetary disk. While L04 also derives an approximate value for the constant $C$, the heuristic nature of the derivation \citep[which is equivalent to a dimensional analysis, as derived by ][]{Johnson06} suggests a better approach would be to fit the amplitude of the perturbations $\gamma$ to the results of MHD simulations.

To compare with the results of G11, we integrated the orbits of an ensemble of 100 planetesimals, started at 5 AUs with zero eccentricity and random initial phase, subject to the stochastic torques described above and no aerodynamic drag and computed the diffusion of their orbital elements. We consider several models, differing only by the turbulent amplitude $\gamma$: model A ($\gamma = 2.5\times 10^{-3}$) represents a fiducial turbulent amplitude appropriate to fully MRI-active disks, while models B1-4 ($\gamma = 2.5\times 10^{-4}$, $10^{-4}$, $5 \times 10^{-5}$, $2.5 \times 10^{-5}$, respectively) have reduced turbulent amplitudes appropriate to the midplane of dead zones.  Since $\alpha$ scales as $\gamma^2$ \citep[where $\alpha$ is the usual viscosity parameter in the Shakura-Sunyaev prescription; ][]{Baruteau10}, the values chosen span two orders of magnitude of turbulent viscosity. Finally, for comparison with M12, we also ran simulations with no turbulence (model C).

 Figure \ref{fig:turb_ecc} shows the growth of the dispersion of the planetesimal eccentricity $\sigma(e)$ as a function of time; as expected from a random-walk process, $\sigma(e) \propto t^{1/2}$. Models B1 and B2 bracket well the growth of the eccentricity dispersion seen by G11 for models with dead zones; we include the reduced values of models B3 and B4 to test the robustness of our results.

We remark that it is likely that the interaction between the planetesimals and the background disk is more complicated than the model presented in this Letter; for instance, we expect the disk to be endowed with some eccentricity \citep[e.g.][]{Marzari08, Marzari12}. However, for the sake of simplicity and to highlight the role of turbulence as an additional factor in dephasing planetesimal orbits, we decided to take the orthogonal approach of ignoring the self-consistent evolution of the disk. In this picture, we assume a circular disk as a maximally accretion-friendly starting point.

\section{Simulations}\label{sec:results}
\subsection{Diffusion of eccentricity and dephasing}

\begin{figure}
\epsscale{1}
\plotone{f2.eps}
\plotone{f3.eps}
\caption{Eccentricity and longitude of pericenter $\varpi - \varpi_B$ as a function of semi-major axis for 10,000 planetesimals distributed between 4 and 8 AU, for runs with no turbulence (black dots), model A (grey dots, active-MRI disk), model B1 (red dots, dead zone) and model B4 (yellow dots, dead zone turbulence amplitude reduced by a factor of 10). }\label{fig:ecclop}
\end{figure}

To visualize the effect of the turbulent fluctuation on a population of planetesimal, we first integrated the trajectories of a swarm of 10,000 planetesimals uniformly distributed between 4 and 8 AUs for four different levels of turbulence: no turbulence, model A (active MRI-disk), model B1 (nominal dead zone turbulence fitting the results of G11) and model B4 (turbulence amplitude reduced by a factor of 10).

Figure \ref{fig:ecclop} shows the eccentricity and longitude of pericenter ($\varpi$) of the planetesimals. Compared to the run without turbulence, both models with levels of turbulence appropriate to a dead zone add a considerable amount of noise on top of the $e$ and $\varpi$-profiles secularly imposed by the central binary. Since the orbits of neighboring planetesimals will not be as collimated, impact velocities will be raised everywhere in the disk. Therefore, we expect that the fraction of accreting impacts will be decreased even when very low levels of turbulence are present.

For the model representing a fully MRI active disk, eccentricities are raised to very high values ($e \approx 5\times 10^{-2}$) and the longitude of pericenter is completely randomized; therefore, the collision velocity can be directly estimated as $\Delta v \approx e v_{kep} > 400$ m/s, comfortably above any velocity threshold resulting in destructive impacts \citep{Stewart09}. Indeed, in this case the turbulent torques completely overwhelm the secular forcing of the central binary, resulting in high impact speeds that are consistent with the results of simulations of planetesimal dynamics in fully turbulent disks around single stars \citep[e.g.][]{Ida08, Nelson10}.

\subsection{Collision statistics}

We subsequently ran full simulations with collision detection between 4 and 10 AU, for each of our models. Given the computational overhead of calculating and updating the turbulent forcing at each time step, we chose to instead integrate the trajectories of a smaller number of planetesimals (1,000) concentrated in annuli centered around 4, 6, 8 and 10 AUs. This approach has the advantage of reducing the running time of our simulations, while simultaneously improving our collision statistics by increasing the impact rate.

\begin{figure}
\epsscale{1}
\plotone{f4.eps}
\caption{Fraction of accreting impacts as a function of semi-major axis (from lighter to darker: models B1-B4 and C.) Each bar represents the percentage of unperturbed impacts; the fraction of perturbed impacts is additionally represented by a solid line.}\label{fig:frac}
\end{figure}

Figure \ref{fig:frac} shows the fraction of accreting impacts as a function of semi-major axis for each of our models. For the sake of completeness, we further classify impacts as ``unperturbed'' ($\Delta v < v_{esc}$, where $v_{esc}$ is the escape velocity) or ``perturbed'' accretion ($v_{esc} < \Delta v < v_{ero}$); for the former, the impact velocity is low enough to allow runaway growth to proceed, therefore allowing the rapid formation of oligarchs. We generously deem a radial bin as a favorable location for planet formation when the fraction of accreting impacts (unperturbed or perturbed) is larger than 50\%.

For models B1 and B2 (with turbulent amplitudes which best fit the results of the 3D MHD simulations of G11), we find that the percentage of accreting impacts decreases dramatically. This is consistent with the noisiness of the pericenter phasing and the high-frequency jitter around the damped value. These short-term oscillations are not damped efficiently by the aerodynamic drag at these large distances, since the drag torque decreases steeply with distance ($\tau_G \propto a^{-3-\beta}$, for planetesimals at the equilibrium forced eccentricity) while the turbulent torque in our model decreases more slowly ($\tau_{T} \propto a^{1-\beta}$), where we take $\Sigma \propto r^{-\beta}$ ($\beta = 1.75$ for the minimum-mass solar nebula model considered in this paper).

For models B3 and B4 (with reduced turbulent amplitudes) we find that, despite the fact that substantial oscillations in eccentricity and longitude of pericenter are still induced, the resulting impact velocities have crossed the accreting threshold. The median impact speed for model B4 is $\approx 6$ m/s, which is in the erosive regime for 1-km planetesimals but allow accretion for 10-km planetesimals.

Finally, in accordance with the simulations of M12, we find that planet formation can proceed undisturbed outside 4 AU when turbulent forcing is switched off.

\section{Discussion}\label{sec:conc}
In this Letter we have investigated planetesimal accretion in the outer parts of a circumbinary disk (4-10 AU). Our simulations indicate that the stochastic forcing of turbulent perturbations will frustrate planetesimal phasing and raise eccentricities, inhibiting planet formation even quite far from the central binary. This result is robust for levels of turbulence that match those observed in realistic MHD simulations. In runs modeling disks fully invaded by MRI turbulence, planetesimal phasing is completely destroyed and planetesimals will collide at speeds that are much higher than fiducial erosive velocities. We then ran simulations fitting the turbulent amplitude corresponding to the stratified disks endowed with a large dead zone obtained in G11, resulting in a reduction by a factor $\approx 10-20$ in the turbulent torque. While the reduced stochastic perturbations would be sufficiently small to allow accretion in a single-star environment, in our circumbinary configuration the fraction of impacts resulting in accreting events was greatly reduced ($\lesssim 1-30\%$ between 4 and 10 AU for models B1 and B2, respectively) due to the increased eccentricity jitter and dephasing. Therefore, it is possible to inhibit planet formation by two different mechanisms: by high eccentricities, differential phasing, and perturbations caused by the self-consistent reaction of the gas disk at small radii, and by stochastic turbulent torques at large radii. This is problematic, since if that picture is accurate, it is not feasible to form the observed planetary census at any realistic distance from the binary. The circumbinary environment is less robust than single stars with regards to dynamical perturbations induced by the disk (regardless of their origins); a further reduced level of turbulence may be required to form planets at all. Indeed, we found that further reduction of the turbulent amplitude by a factor of 2-5 was necessary to make planet formation viable again in the outer disk.

Several uncertainties in our model may still provide some room for allowing planet formation despite the perturbing effects of turbulence. Our analytic prescription, while attempting to match MHD simulations, only provides a ``0-th'' order description of gravitational perturbations induced by torques. For instance, we remark that we tuned the amplitude of turbulent torques based on the evolution of swarm of planetesimals in a single MHD simulation, and rescaled the turbulent amplitude at each radial location according to the numerical prescription of Section \ref{sec:setup}; however, the output of these models will depend on the assumed magnetic field, the ionization level and possibly the resolution of the shearing box.
The detailed radial dependence and even extent of the dead zone remain uncertain and depend on the assumed model \citep[e.g.][; if the dead zone extends to only a few AUs, then fully active turbulence may still play a role in the outer disk and even shut off planetesimal accretion completely]{Matsumura05, Terquem08, Flaig12}. Effectively, using the torque scaling of Equation \ref{eqn:torque} implicitly assumes that the turbulent torque is well described by the fully turbulent model of L04, attenuated by a constant factor at each radial location. Only self-consistent (including the time-dependent potential of the binary), computationally expensive global simulations could, in principle, inform the model presented in this Letter. However, since we have considered various values of $\gamma$ for each independent annulus, we are reasonably covering a number of possibilities for the radial dependence of the turbulence amplitude. Our results at a given radius should only be considered suggestive of further challenges to planet formation.

Finally, we remark that, similarly to M12, we have chosen to study planetesimal growth through mutual collisions in the 1-10 km planetesimal regime. A more sophisticated approach following the evolution of the size distribution of the planetesimals (adding a significant computational overhead) may be warranted; in the latter approach, coevolution with a disk of small dust might make accretion possible closer to the binary \citep[assuming a high dust accretion efficiency; ][]{Paardekooper12}. Even in this setup, some amount of migration is needed to bring the fully formed core (or embryos) to the current location. The assumption of a static gas background is not realistic, especially in the inner disk; in that case, self-consistent gas dynamics will act to further increase impact velocities in the inner disk \citep{Paardekooper08}. We plan to incorporate the full hydrodynamical evolution of the gas disk in a follow-up paper using the SPH module of our code.

Our simulations indicate that planetesimal accretion in the 1-10 km range will be inhibited everywhere in the disk, and could only proceed if the strength of the turbulent torques is reduced from our fiducial value, or the initial planetesimal population is comprised of bigger objects. The latter scenario is particularly appealing, as it could allow planet formation to proceed in other highly dynamically disturbed environments \citep[e.g.][]{Thebault11}. Indeed, our fiducial model for dead zone turbulence becomes accretion-friendly for a planetesimal size spectrum spanning 10-100 km size. Recent simulations of planetesimal formation in weakly turbulent disks show that massive bound clumps are formed rapidly from meter-sized boulders within pressure bumps, for typical nebula parameters \citep[e.g. ][]{Johansen07, Johansen11}. Such bound clumps will likely result in planetesimals with sizes comparable to at least a substantial fraction of the dwarf planet Ceres. We suggest that a primordial population of large planetesimals might be crucial to proceed with planet formation in highly perturbed environments. However, further simulations will be required to assess whether the formation of large clumps is robust to the dynamical perturbations of a binary companion.

\acknowledgments
S.M. acknowledges insightful discussions with G. Laughlin and N. Haghighipour, a useful critique from the anonymous referee, as well as support from the W. J. McDonald Postdoctoral Fellowship.

\bibliographystyle{apj}


\section{Introduction}
In this note we recall some classical results on the structure of the subdifferential of first order one-homogeneous functionals, and
we give new regularity results which extend and precise previous work by G. Anzellotti \cite{Anzelotti,anzeltrace,anzeleuler}.

Given an open set $\Omega\subset\mathbb{R}^d$  with Lipschitz boundary, and a function $u\in C^1(\Omega)\cap BV(\Omega)$, we consider the functional
\[
J(u)\ :=\ \int_\Omega F(x,Du)
\]
where $F:\Omega\times\mathbb{R}^d\to [0,+\infty)$ is continuous,
and $F(x,\cdot)$ is a smooth and uniformly convex positively
one-homogeneous functional on $\mathbb{R}^d$
for all $x\in \Omega$. The functional $J$ can be canonically
relaxed to the whole of $BV(\Omega)$ (see \cite[Section 5.5]{AFP})
and we still write, in analogy with the notation commonly used for the total variation,
$J(u)=\int_\Omega F(x,Du)$ for $u\in BV(\Omega)$,
where $F(x,Du)$ is in general a Radon measure in $\Omega$ defined by $F(x,Du):=F\left(x,\frac{Du}{|Du|}\right)|Du|$ with $\frac{Du}{|Du|}$ the Radon-Nikodym derivative of $Du$ with respect to $|Du|$.

Since $BV(\Omega)\subset L^{d/(d-1)}(\Omega)$, it is
natural to consider $J$ as a convex, l.s.c.~function on the whole of
$L^{d/(d-1)}(\Omega)$,
with value $+\infty$ when $u\not\in BV(\Omega)$.
In this framework,
for any $u\in L^{d/(d-1)}(\Omega)$ such that $J(u)<+\infty$,
that is $u\in BV(\Omega)$,
we can define the subdifferential of $J$
at $u$, in the duality $(L^{d/(d-1)},L^d)$, as
\begin{equation*
\partial J(u)\ :=\ \left\{
g\in L^d(\Omega)\,:\, J(v)\ge J(u)\,+\,\int_\Omega g(v-u)\,dx\ \forall
v\in L^{d/(d-1)}(\Omega) \right\}.
\end{equation*}
Notice that a function $g\in L^d(\Omega)$ belongs to $\partial J(u)$
if and only if $u$ is a minimizer of the functional $J(v)-\int_\Omega gv dx$
among all $v\in L^{d/(d-1)}(\Omega)$.

The goal of this paper is to investigate the particular structure
of the functions $u$ and $g$, when $g\in \partial J(u)$.
Let
\[
F^*(x,z) := \sup_{w\in\mathbb{R}^d} z \cdot w - F(x,w)
\]
be the Legendre-Fenchel or convex conjugate of $F$.
Notice that $F^*$ takes values in $\{0,+\infty\}$,
and $F^*(x,z)=0$ if and only if $F^\circ(x,z)\le 1$, where
$F^\circ$ denotes the convex polar of $F$ defined in~\eqref{polar} below.

Given $u\in L^{d/(d-1)}$, the functional $J(u)$ can also be expressed by duality as
\[
J(u)\ =\ \sup\left\{
-\int_\Omega u \textup{div}\, z \,dx\ :\ z\in C_c^\infty(\Omega;\mathbb{R}^d)\,,\ F^*(x,z(x))=0
\ \forall x\in\Omega
\right\}.
\]
It follows that a function $g\in\partial J(u)$ has necessarily the form
$g=-\textup{div}\, z$, for some vector field
$z\in L^\infty(\Omega;\mathbb{R}^d)$ with $F^*(x,z(x))=0$ a.e.~in $\Omega$.
Since by a formal integration by parts one gets $$z\cdot Du= F(x,Du),$$
two natural questions arise:
\begin{itemize}
\item in what sense is this relation true?
\item can one assign a precise value to $z$ on the support of the measure $Du$?
\end{itemize}

The first question has been answered by
Anzelotti in the series of papers~\cite{Anzelotti,anzeltrace,anzeleuler}.
However, for the particular vector fields we are interested in,
we can be more precise and obtain pointwise properties of $z$
on the level sets of the function $u$. Indeed, we shall show
that $z$ has a pointwise meaning on all level sets of
$u$, up to $\mathcal{H}^{d-1}$-negligible sets (which can be
much more than $|Du|$-a.e.,
as illustrated by the function $u =\sum_{n=1}^{+\infty} 2^{-n}\chi_{(0,x_n)}$,
defined in the interval $(0,1)$, with $(x_n)$ a dense sequence in that
interval).

We will therefore focus on the properties of the vector fields
$z\in L^\infty(\Omega,\mathbb{R}^d)$ such that $F^*(x,z(x))=0$ a.e. in $\Omega$
and $g=-\textup{div}\, z \in  L^d(\Omega)$, and such that there exists a function
$u$ such that for any $\phi\in C_c^\infty(\Omega)$,
\[
-\int_\Omega\textup{div}\, z \,u\phi\,dx\ =\ \int_\Omega u\,z\cdot\nabla\phi\,dx
+\int_\Omega \phi F(x,Du)\,.
\]
In particular, one checks easily that $u$  minimizes the functional
\begin{equation}\label{minsubgrad}
\int_\Omega F(x,Du) - \int_\Omega g u \,dx
\end{equation}
among perturbations with compact support in $\Omega$.
Conversely, given $g\in L^d(\Omega)$ with $\|g\|_{L^d}$ sufficiently small,
there exist functions $u$ which
minimize~\eqref{minsubgrad} under various types of boundary conditions,
and corresponding fields $z$.

This kind of functionals appears in many contexts including image processing and plasticity \cite{Andreu,Temam}.
Notice also that, by the Coarea Formula \cite{AFP}, it holds
\begin{equation*
\int_\Omega F(x,Du) - \int_\Omega g u\,dx
= \int_\mathbb{R} \left( \int_{\partial^* \{u>s\}} F(x,\nu) \,d\mathcal{H}^{d-1}(x)
- \int_{\{u>s\}} g\,dx\right)\,ds\,,
\end{equation*}
where $\nu$ is the unit normal to $\{u>s\}$,
and one can show (see for instance~\cite{CMCF}) that
any level set of the form
$\{ u>s\}$ or $\{ u\ge s\}$ is a minimizer of the geometric functional
\begin{equation}\label{geomfun}
E\mapsto \int_{\partial^* E} F(x,\nu)\,d\mathcal{H}^{d-1}(x) - \int_E g\,dx\,.
\end{equation}
defined for sets $E$ of finite perimeter.
The canonical example of such functionals is given by the total variation,
corresponding to $F(x,Du)=|Du|$.
In this case, \eqref{geomfun} boils down to
\begin{equation}\label{geomper}
P(E)-\int_E g \,dx.
\end{equation}

In \cite{BarGonTa}, it is shown that every set with finite perimeter in $\Omega$
is a minimizer of \eqref{geomper} for some $g\in L^1(\Omega)$.
If $F$ is even in $\nu$ and when $g\in L^d(\Omega)$, the boundary $\partial E$ is only of class $C^{0,\alpha}$ out of a singular set (see \cite{AmbPao}). However, if $g\in L^p(\Omega)$ with $p>d$,
and $E$ is a minimizer of \eqref{geomfun},
then $\partial E$ is locally $C^{1,\alpha}$ for some $\alpha>0$, out of a closed singular set of zero $\mathcal{H}^{d-3}$-measure \cite{ASS,SchoenSimon} (some
regularity assumption on $F$ is required, see also Remark~\ref{remDu} below).
 Since the Euler-Lagrange equation of \eqref{geomfun} relates $z$ to the normal of $E$,
understanding the regularity of $z$ is closely related to understanding the regularity of $\partial E$.

\smallskip

Our main result is that the Lebesgue points of $z$ correspond to regular points of $\partial \{u>s\}$ or $\partial \{u\ge s\}$ (Theorem \ref{zfixPL}),
and that the converse is true in dimension $d\le 3$ (Theorem \ref{propLebesgue}).

\section{Preliminaries}

\subsection{BV functions}
We briefly recall the definition of function of bounded variation and set of finite perimeter. For a complete presentation we refer to \cite{AFP}.

\begin{definition}
Let $\Omega$ be an open set of $\mathbb{R}^d$, we say that a function $u\in L^1(\Omega)$ is a function of bounded variation if
\[\int_\Omega |Du|:= \sup_{\stackrel{z \in \mathcal{C}^1_c(\Omega)}{|z|_\infty \leq 1}} \int_{\Omega} u \ \textup{div}\, z \ dx< +\infty.\]
We denote by $BV(\Omega)$ the set of functions of bounded variation in $\Omega$ (when $\Omega=\mathbb{R}^d$ we simply write $BV$ instead of $BV(\mathbb{R}^d)$).
We say that a set $E\subset \mathbb{R}^d$ is of finite perimeter if its characteristic function $\chi_E$ is of bounded variation and denote its perimeter in an open set $\Omega$ by $P(E,\Omega):=\int_{\Omega} |D \chi_E|$, and write simply $P(E)$ when $\Omega=\mathbb{R}^d$.
\end{definition}

\begin{definition}
Let $E$ be a set of finite perimeter and let $t\in [0;1]$. We define
\[
E^{(t)}:= \left\{ x \in \mathbb{R}^d \,:\, \lim_{r\downarrow 0} \frac{|E\cap B_r(x)|}{|B_r(x)|}=t\right\}.
\]
We denote by $\partial E:= \left(E^{(0)}\cup E^{(1)}\right)^c$ the measure theoretical boundary of $E$.
Let $\textup{spt}\,(|D \chi_E|)$ be the support of the measure $|D \chi_E|$:
we define the reduced boundary of $E$ by:
\[
\partial^* E:= \left\{ x \in \textup{spt}\,(|D \chi_E|) \,: \, \nu^E(x):= \lim_{r \downarrow 0} \frac{ D \chi_E (B_r(x))}{|D \chi_E|(B_r(x))} \; \textrm{exists and } \; |\nu^E(x)|=1 \right\}\subset E^{\left(\frac 1 2\right)}.
\]
The vector $\nu^E(x)$ is the measure theoretical inward normal to the set $E$.
\end{definition}
\begin{proposition}
If $E$ is a set of finite perimeter then $D\chi_E=\nu^E\  \mathcal{H}^{d-1} {\mathchoice {\>\LM7\>}{\>\LM7\>}{\,\LM5\,}{\,\LM{3.35}\,}}\partial^* E$, $P(E)=\mathcal{H}^{d-1}(\partial^* E)$ and $\mathcal{H}^{d-1}(\partial E \setminus \partial^* E)=0$.
\end{proposition}

\begin{definition}
 We say that $x$ is an approximate jump point of $u\in BV(\Omega)$
if there exist  $\xi \in \mathbb{S}^{d-1}$
and distinct $a, b \in \mathbb{R}$
 such that
\[ \lim_{\rho \to 0} \frac{1}{|B_\rho^+(x,\xi)|} \int_{B_\rho^+(x,\xi)} |u(y)-a| \ dy =0 \quad \textrm{ and } \quad  \lim_{\rho \to 0} \frac{1}{|B_\rho^-(x,\xi)|} \int_{B_\rho^-(x,\xi)} |u(y)-b| \ dy =0,\]
where $B_\rho^\pm(x,\xi):= \{ y \in B_\rho(x) : \pm (y-x)\cdot{\xi}>0 \}.$ Up to a permutation of $a$ and $b$ and a change of sign of $\xi$, this characterize the triplet $(a,b,\xi)$ which is then denoted by $(u^+,u^-, \nu_u)$. The set of approximate jump points is denoted by $J_u$.
\end{definition}

The following proposition can be found in \cite[Proposition 3.92]{AFP}.
\begin{proposition}\label{Ju}
 Let $u\in BV(\Omega)$. Then, defining
\[\Theta_u:=\{ x\in \Omega \, : \, \liminf_{\rho \to 0} \rho^{1-d} |Du|(B_\rho(x))>0\},\]
there holds $J_u\subset \Theta_u$ and $\mathcal{H}^{d-1}(\Theta_u \setminus J_u)=0$.
\end{proposition}

\subsection{Anisotropies}
Let $F(x,p):\Omega \times \mathbb{R}^d\to [0,+\infty)$ be a continuous functions, which is
convex and positively one-homogeneous in the second variable:
$F(x,\lambda p)=\lambda F(x,p)$ for all $\lambda>0,x,p$;
and such that there exists $c_0>0$ with
\[c_0 |p|\le F(x,p)\le \frac{1}{c_0} |p| \qquad \forall (x,p)\in \mathbb{R}^d\times \mathbb{R}^d.\]
We say that $F$ is uniformly elliptic if for some $\delta>0$, the function
$p\mapsto F(p)-\delta |p|$ is still a convex function.  We define the polar function of $F$ by
\begin{equation}\label{polar}
F^\circ(x,z):=\sup_{\{F(x,p)\le1\}} z\cdot p
\end{equation}
so that $(F^\circ)^\circ=F$.  It is easy to check that
$[F(x,\cdot)^2/2]^*= F^\circ(x,\cdot)^2/2$, where as before the $^*$ denotes
the convex conjugate with respect to the second variable.
In particular, if differentiable,
$F(x,\cdot)\nabla_p F(x,\cdot)$ and $F^\circ(x,\cdot)\nabla_z F^\circ(x,\cdot)$
are inverse monotone operators.
Also, one has that $F^*(x,z)=0$ if and only if $F^\circ(x,z) \le 1$,
and $F^*(x,z)=+\infty$ else.

If $F(x,\cdot)$ is differentiable then, for every $p\in \mathbb{R}^d$,
\[
F(x,p)=p\cdot \nabla_p F(x,p) \qquad {\rm (Euler's\ identity)}
\]
and
\[ z\in \{F^\circ(x,\cdot)\le 1\} \textrm{ with } p\cdot z=F(x,p)\ \Longleftrightarrow \  z=\nabla_p F(x,p).\]
If $F$ is elliptic and of class $\mathcal{C}^2(\mathbb{R}^d \times \mathbb{R}^d\setminus\{0\})$, then $F^\circ$ is also elliptic and $\mathcal{C}^2(\mathbb{R}^d \times \mathbb{R}^d\setminus\{0\})$. We will then say that $F$ is a smooth elliptic anisotropy.
Observe that, in this case, the function $F^2/2$ is also uniformly $\delta^2$-convex
(this follows from the inequalities $D^2F(x,p)\ge \delta/|p|(I-p\otimes p/|p|^2)$
and $F(x,p)\ge \delta |p|$). In particular, for every $x,y,z \in \mathbb{R}^d$, there holds
\begin{equation}\label{strictconvF}F^2(x,y)-F^2(x,z)\ge 2\left(F(x,z)\nabla_p F(x,z)\right)\cdot (y-z) + \delta^2|y-z|^2,
\end{equation}
and a similar inequality holds for $F^\circ$.
We refer to \cite{schneider} for general results on convex norms and convex bodies.

\subsection{Pairings between measures and bounded functions}
Following \cite{Anzelotti} we define a generalized trace $[z,Du]$
for functions $u$ with bounded variation and bounded vector fields $z$ with divergence in $L^d$.
\begin{definition}[Anzelotti~\cite{Anzelotti}]
Let $\Omega$ be an open set with Lipschitz boundary,  $u\in BV(\Omega)$ and
 $z \in L^{\infty}(\Omega,\mathbb{R}^d)$ with $\textup{div}\, z\in L^d(\Omega)$.
We define the distribution $[z,Du]$ by
 \[ \langle [z,Du], \psi \rangle= -\int_\Omega u \,\psi\, \textup{div}\, z \, dx - \int_\Omega u \,  z \cdot \nabla \psi  \, dx
 \qquad \forall \psi \in \mathcal{C}^\infty_c(\Omega).
 \]
\end{definition}

\begin{proposition}[Anzelotti~\cite{Anzelotti}]
The distribution $[z,Du]$ is a bounded Radon measure on $\Omega$ and if $\nu$ is the inward unit normal to $\Omega$, there exists a function $[z,\nu]\in L^\infty(\partial \Omega)$ such that the generalized Green's formula holds,
 \[ \int_\Omega [z,Du]=- \int_\Omega u \textup{div}\, z \, dx -\int_{\partial \Omega} [z, \nu] u\,d \mathcal{H}^{d-1}.\]
 The function $[z,\nu]$ is the generalized (inward) normal trace of $z$ on $\partial \Omega$.
\end{proposition}

Given $z \in L^{\infty}(\Omega;\mathbb{R}^d)$, with $\textup{div}\, z\in L^d(\Omega)$, we can
also define the generalized trace of $z$ on $\partial E$, where $E$ is a set of locally finite perimeter.
Indeed, for every bounded open set $\Omega$ with Lipschitz boundary, we can define as above the measure $[z,D\chi_E]$ on $\Omega$.
Since this measure is absolutely continuous with respect to $|D\chi_E|= \mathcal{H}^{d-1} {\LL} \partial^* E$ we have
\[
[z,D\chi_E]= \psi_z(x) \mathcal{H}^{d-1}{\LL} \partial^* E
\]
with $\psi_z \in L^{\infty}(\partial^* E)$ independent of $\Omega$.
We denote by $[z,\nu^E]:=\psi_z$ the generalized (inward) normal trace of $z$ on $\partial E$.
If $E$ is a bounded set of finite perimeter, by taking $\Omega$ strictly containing $E$,
we have the generalized Gauss-Green Formula
\[\int_E \textup{div}\, z \, dx =-\int_{\partial^* E} [z,\nu^E] d \mathcal{H}^{d-1}.\]
Anzellotti proved the following alternative definition of $[z,\nu^E]$ \cite{anzeltrace,anzeleuler}

\begin{proposition}[Anzelotti~\cite{anzeltrace,anzeleuler}]
\label{convcylinder}
Let $(x,\alpha)\in \mathbb{R}^d\times \mathbb{R}^d\setminus\{0\}$.
For any $r>0$, $\rho>0$ we let
$$
C_{r,\rho}(x,\alpha):=\{ \xi \in \mathbb{R}^d \, : \,
\left|(\xi-x)\cdot \alpha\right| <r,
\left|(\xi-x)-[(\xi-x)\cdot \alpha] \alpha\right|<\rho\}.
$$
There holds
\[ [z,\alpha](x) =\lim_{\rho\to 0}\lim_{r\to 0} \  \frac{1}{2r \omega_{d-1} \rho^{d-1}} \int_{C_{r,\rho}(x,\alpha)} z \cdot \alpha\, dy \]
where $\omega_{d-1}$ is the volume of the unit ball in $\mathbb{R}^{d-1}$.
\end{proposition}

\section{The subdifferential of anisotropic total variations}

\subsection{Characterization of the subdifferential}

The  following characterization of the  subdifferential of $J$ is classical and readily follows for example from the representation
formula~\cite[(4.19)]{DMBou}.
\begin{proposition}\label{defcalibration}
Let $F$ be a smooth elliptic anisotropy and $g\in L^d(\Omega)$ then $u$ is a local minimizer of \eqref{minsubgrad} if and only if there exists $z\in L^\infty(\Omega)$ with $\textup{div}\, z =g$, $F^*(x,z(x))=0 $ a.e. and
\[ [z,Du]= F(x,Du),\]
where the equality holds in the sense of measures.
Moreover, for every $t\in \mathbb{R}$, for the set $E=\{u>t\}$ there holds $[z,\nu^E]=F(x,\nu^E) $ $\mathcal{H}^{d-1}$-a.e. on $\partial E$. We will say that such a vector field is a calibration of the set $E$ for the minimum problem~\eqref{geomfun}.
\end{proposition}

\begin{remark}\rm
In \cite{Anzelotti}, it is proven that if $z_\rho(x):=\frac{1}{|B_\rho(x)|}\int_{B_\rho(x)} z(y) \ dy$, then $z_\rho\cdot\nu^E$  weakly* converges to $[z, \nu^E]$ in $L^\infty_{loc}(\mathcal{H}^{d-1}{\mathchoice {\>\LM7\>}{\>\LM7\>}{\,\LM5\,}{\,\LM{3.35}\,}}\partial^* E)$. Using \eqref{strictconvF} it is then possible to prove that if $z$ calibrates $E$ then $z_\rho$ converges to $\nabla_p F(x,\nu^E)$ in $L^2(\mathcal{H}^{d-1}{\mathchoice {\>\LM7\>}{\>\LM7\>}{\,\LM5\,}{\,\LM{3.35}\,}} \partial^*E)$ yielding that up to a subsequence,
$z_{\phi(\rho)}$ converges $\mathcal{H}^{d-1}$-a.e. to $\nabla_p F(x, \nu^E)$. Unfortunately this is still a very weak statement since it is a priori impossible to recover from this the convergence of the full sequence $z_\rho$.
\end{remark}

The main question we want to investigate now is whether we can give a classical meaning to $[z,\nu^E]$ (that is understand if $[z,\nu^E]=z\cdot \nu^E$). We observe that a priori the value of $z$ is not well defined on $\partial E$ which has zero Lebesgue measure.

We let $S:=\textup{spt}\,(Du)\subset\Omega$ be the support of the measure
$Du$, that is, the smallest closed set
in $\Omega$ such that $|Du|(\Omega\setminus S)=0$. We will show that essentially
in $S$, $z$ is well-defined, as soon as $g\in L^d(\Omega)$.

The next result is classical, for a proof we refer to \cite{Massari,GoMa}.
\begin{lemma}[Density estimate]\label{lemdens}
There exists $\rho_0>0$ (depending on $g\in L^d(\Omega)$)
and a constant $\gamma>0$ (which depends only on $d$), such that for
any $B_\rho(x)\subset \Omega$
with $\rho\le\rho_0$, and any level set $E$ of $u$ (that is,
$E\in \left\{\{u>s\},\{u\ge s\},\{u<s\},\{u\le s\},\, s\in \mathbb{R}\right\}$), if
$|B_\rho(x)\cap E|<\gamma |B_\rho(x)|$
then $|B_{\rho/2}(x)\cap E|=0$.
As a consequence, $E^{0}$ and $E^{1}$ are open, $\partial E$ is the topological
boundary of $E^{1}$,
and (possibly changing slightly $\gamma$)
if $x\in \partial E$, then $\mathcal{H}^{d-1}(\partial E\cap B_\rho(x))
\ge \gamma \rho^{d-1}$.
\end{lemma}
This result is not true anymore if $g\not\in L^d(\Omega)$~\cite{GoMa}.
If $\partial \Omega$ is Lipschitz, it is true up to the boundary.
\medskip

\begin{corollary}
It follows that $u\in L^\infty_{loc}(\Omega)$\, and $u\in C(\Omega\setminus \Theta_u)$.
\end{corollary}
\begin{proof}
 For any ball $B_\rho(x)\subset \Omega$ and
 $\inf_{B_{\rho/2}(x)}u<a<b<\sup_{B_{\rho/2}(x)}u$, one has
\[
+\infty\,>\, |Du|(B_\rho(x))\,\ge \,\int_{a}^b P(\{u>s\},B_\rho(x))\,ds
\,\ge\, (b-a)\gamma\left(\frac{\rho}{2}\right)^{d-1}\,,
\]
so that $osc_{B_{\rho/2}(x)}(u)$ must be bounded and thus $u\in L^\infty_{loc}(\Omega)$. Moreover, if $x \in \Omega\setminus \Theta_u$ we find that
 $\lim_{\rho \to 0} osc_{B_{\rho}(x)}(u)=0$ so that $u$ is continuous at the point $x$.
\end{proof}

We remark that if sets $(E_n)_n$ satisfy the density estimate of Lemma~\ref{lemdens}
and converge in $L^1$ to some limit set, then one easily deduces that
the convergence also holds in the Hausdorff (or Kuratowski, is the sets
are unbounded) sense. Applying this principle to the level sets of $u$,
we find that all points in the support
of $Du$ must be on the boundary of a  level set of $u$:
\begin{proposition}\label{proplevelset}
For any $x\in S$, there exists $s\in \mathbb{R}$ such that either $x\in \partial\{u >s\}$
or $x\in \partial \{u\ge s\}$.
\end{proposition}
\begin{proof} First, if $x\not\in S$ then $|Du|(B_\rho(x))=0$
for some $\rho>0$ and clearly $x$ cannot be on the boundary of a level
set of $u$. On the other hand, if $x\in S$, then for any ball $B_{1/n}(x)$
($n$ large) there is a level $s_n$ (uniformly bounded) with
$\partial \{u>s_n\}\cap B_{1/n}(x)\neq \emptyset$ and by Hausdorff convergence,
we deduce that either $x\in\partial\{u>s\}$ or  $x\in\partial \{u\ge s\}$
where $s$ is the limit of the sequence $(s_n)_n$ (which must actually converge).
\end{proof}
\medskip
The following stability property is classical (see e.g. \cite{CGN}).
\begin{proposition}\label{stabilitycalib}
Let $E_n$ be
local minimizers of~\eqref{geomfun}, with a function $g=g_n\in L^d(\Omega)$,
and converging in the $L^1$-topology to a set $E$.
Assume that the sets $E_n$ are calibrated by $z_n$ (see Prop.~\ref{defcalibration}),
that $z_n\stackrel{*}{\rightharpoonup} z$ weakly-$*$ in $L^\infty$
and $g_n\to g=-\textup{div}\, z\in L^d(\Omega)$, in $L^1(\Omega)$ as $n\to\infty$.
Then $z$ calibrates $E$, which is
thus also a minimizer of~\eqref{geomfun}.
\end{proposition}
Let us observe that, if $z_n\stackrel{*}{\rightharpoonup}z$
and $F^\circ(x,z_n(x))\le 1$ a.e. in $\Omega$, then in the limit one still gets
$F^\circ(x,z(x))\le 1$ a.e. in $\Omega$, thanks to the continuity of $F$
and the convexity in the second variable.

\subsection{The Lebesgue points of the calibration.}

The next result shows that the regularity of the calibration $z$
implies some regularity of the calibrated set.

\begin{theorem}\label{zfixPL}
Let $E=\{u>t\}$ or $E=\{u\ge t\}$, and
let $\bar x\in \partial E$ be a Lebesgue point of $z$.
Then, $\bar x  \in \partial^*E$ and
\begin{equation}\label{zeqnu}
z(\bar x)=\nabla_p F(\bar x,\nu^E(\bar x)).
\end{equation}
\end{theorem}

\begin{proof}
We reason as in \cite[Th. 4.5]{CGN} and let $z_\rho(y) := z(\bar x+\rho y)$.
Since $\bar x$ is a Lebesgue point of $z$, we have
that $z_\rho\to\bar z$ in $L^1(B_R)$, hence
also weakly-$*$ in $L^\infty(B_R)$ for any $R>0$,
where $\bar z\in \mathbb{R}^d$ is a constant vector.

We let $E_\rho := (E-\bar x)/\rho$ and  $g_\rho(y):=g(\bar x + \rho y)$ (so that $\textup{div}\, z_\rho=\rho g_\rho$).
Observe that $E_\rho$ minimizes
\[
\int_{\partial^* E_\rho\cap B_R} F(\bar x+\rho y,\nu^{E_\rho}(y))\,d\mathcal{H}^{d-1}(y)
+ \rho\int_{E_\rho\cap B_R} g_\rho(y)\,dy\,,
\]
with respect to compactly supported
perturbations of the set (in the fixed ball $B_R$).
Also,
\[
\|\rho g_\rho\|_{L^d(B_R)}\ = \ \|g\|_{L^d(B_{\rho R})}
\ \stackrel{\rho\to 0}{\longrightarrow}\ 0\,.
\]
By Lemma~\ref{lemdens}, the sets $E_\rho$ (and the boundaries $\partial E_\rho$)
satisfy uniform density bounds, and hence are compact with respect
to both local $L^1$ and Hausdorff convergence.

Hence, up to extracting a subsequence, we can assume that $E_\rho\to \bar E$, with $0\in \partial \bar E$.
Proposition~\ref{stabilitycalib} shows that $\bar z$ is a calibration
for the energy $\int_{\partial \bar E\cap B_R} F(\bar x, \nu^{\bar E}(y))\,d\mathcal{H}^{d-1}(y)$,
and that $\bar E$ is a minimizer calibrated by $\bar z$.

It follows that $[\bar z,\nu^{\bar E}]=F(\bar x,\nu^{\bar E}(y))$ for
$\mathcal{H}^{d-1}$-a.e.~$y$ in $\partial \bar E$, but since $\bar z$ is a constant,
we deduce that $\bar E=\{ y\cdot \bar\nu\ge 0\}$ with
$\bar \nu/F(\bar x,\bar \nu)=\nabla_p F^\circ(\bar x,\bar z)$\footnote{
We use here that $F(\bar x,\cdot)\nabla F(\bar x,\cdot)=
[F^\circ(\bar x,\cdot)\nabla F^\circ(\bar x,\cdot)]^{-1}$, so that
$\bar z=\nabla F(\bar x,\nu^{\bar E}(y))$ implies both
$F^\circ(\bar x,\bar z)=1$ and
$\nu^{\bar E}(y)/F(\bar x,\nu^{\bar E})(y)=\nabla F^\circ(\bar x,\bar z)$}.
In particular the limit $\bar E$
is unique, hence we obtain the global convergence of $E_\rho\to \bar E$, without passing to a subsequence.

We want to deduce that $\bar x\in\partial^* E$, with
$\nu^E(\bar x)=F(\bar x,\nu^{E}(\bar x))\nabla_p F^\circ(\bar x,\bar z)$, which is equivalent to \eqref{zeqnu}.
The last identity is obvious from
the arguments above, so that we only need to show that
\begin{equation}\label{reducbound}
\lim_{\rho\to 0} \frac{ D\chi_{E_\rho}(B_1) }{ |D\chi_{E_\rho}|(B_1) }\ =\ \bar \nu\,.
\end{equation}
Assume we can show that
\begin{equation}\label{convenergy}
\lim_{\rho\to 0} |D\chi_{E_\rho}|(B_R)\ =\ |D\chi_{\bar E}|(B_R)
\ \left(\ =\ \omega_{d-1}R^{d-1}\right)
\end{equation}
for any $R>0$, then for any $\psi\in C_c^\infty(B_R;\mathbb{R}^d)$ we would get
\begin{multline*}
\frac{1}{|D\chi_{E_\rho}|(B_R)}\int_{B_R} \psi\cdot D\chi_{E_\rho}
\ =\ -\frac{1}{|D\chi_{E_\rho}|(B_R)}\int_{B_R\cap E_\rho} \textup{div}\, \psi(x)\,dx
\\ \longrightarrow\ -\frac{1}{|D\chi_{\bar E}|(B_R)}\int_{B_R\cap \bar E} \textup{div}\, \psi(x)\,dx
\ =\ \frac{1}{|D\chi_{\bar E}|(B_R)}\int_{B_R} \psi\cdot D\chi_{\bar E}
\end{multline*}
and deduce that the
measure $D\chi_{E_\rho}/(|D\chi_{E_\rho}|(B_R))$  weakly-$*$ converges to
$D\chi_{\bar E}/(|D\chi_{\bar E}|(B_R))$.
Using again \eqref{convenergy}), we then obtain that
\begin{equation}\label{aeR}
\lim_{\rho\to 0} \frac{ D\chi_{E_\rho}(B_R) }{ |D\chi_{E_\rho}|(B_R) }\ =\ \bar \nu
\end{equation}
for almost every $R>0$.
Since $D\chi_{E_\rho}(B_{\mu R})/(|D\chi_{E_\rho}|(B_{\mu R}))=D\chi_{E_{\rho/\mu}}(B_{R})/(|D\chi_{E_{\rho/\mu}}|(B_{R}))$ for any $\mu>0$,
\eqref{aeR} holds in fact for any $R>0$ and \eqref{reducbound} follows, so that
$\bar x\in \partial^* E$.

It remains to show~\eqref{convenergy}. First, we observe that,
by minimality of $E_\rho$ and $\bar E$ plus the Hausdorff convergence of $\partial E_\rho$ in balls,
we can easily show the convergence of the energies
\begin{multline*}
\lim_{\rho\to 0} \int_{\partial E_\rho\cap B_R}  F(\bar x+\rho y,\nu^{E_\rho}(y))\,d\mathcal{H}^{d-1}(y)
+ \rho\int_{E_\rho\cap B_R} g_\rho(y)\,dy
\\ =\ \int_{\partial \bar E\cap B_R}
F(\bar x,\nu^{\bar E}(y))\,d\mathcal{H}^{d-1}(y)
\end{multline*}
and, by the continuity of $F$,
\begin{equation}\label{convenergyrho}
\lim_{\rho\to 0} \int_{\partial E_\rho\cap B_R}  F(\bar x,\nu^{E_\rho}(y))\,d\mathcal{H}^{d-1}(y)
\ =\ \int_{\partial \bar E\cap B_R}
F(\bar x,\nu^{\bar E}(y))\,d\mathcal{H}^{d-1}(y)\,.
\end{equation}
Then, \eqref{reducbound} follows from Reshetnyak's continuity
theorem where, instead of using the Euclidean norm as reference norm, we use the uniformly convex function $F(\bar x, \cdot)$ and the convergence of the measures $F(\bar x, D\chi_{E_\rho})$ to $F(\bar x,D \chi_{\bar E})$ (see  \cite{Resh,CGN}).

\end{proof}

\begin{corollary}
For any $x\in S$ let $E^x \in \left\{\{u>u(x)\},\{u\ge u(x)\}\right\}$
be the upper level set of $u$ such
that $x\in\partial E^x$. Then, the equality
\begin{equation}\label{euleru}
 z(x)=\nabla_p F\left(x,\frac{D\chi_{E^x}}{|D\chi_{E^x}|}(x)\right)
\end{equation}
holds Lebesgue a.e. in $S=\textup{spt}\,(Du)$.
\end{corollary}

\begin{remark}\rm
In the inhomogeneous isotropic case $F(x,p)=a(x)|p|$, with $a(\cdot)$ periodic,
a similar result has been proved by
Auer and Bangert in \cite[Th. 4.2]{AB}. As a consequence they obtain
differentiability properties of the so-called stable norm associated to
the functional $J$ (see also \cite{CGN} for the anisotropic
version of their result).
\end{remark}

In dimension $2$ and $3$ we can also show the reverse implication, proving that regular points of the boundary correspond to Lebesgue points of the calibration. The idea is to show that the parameters $r,\,\rho$ in Proposition \ref{convcylinder} can be taken of the same order.

\begin{theorem} \label{propLebesgue}
Assume the dimension is $d=2$ or $d=3$.
Let $x,s$ be as  in Proposition~\ref{proplevelset}, $E$ be a minimizer of \eqref{geomfun} and assume
$x\in \partial^*E$. Then $x$ is a Lebesgue point of $z$
and \eqref{zeqnu} holds at $x$.
\end{theorem}

\begin{proof}
We divide the proof into two steps.

{\it Step 1.} We first consider anisotropies $F$ which are not depending on the $x$ variable.  Without loss of generality we assume $x=0$.
By assumption, there exists the limit
\begin{equation}
\ov{\nu} \ :=\ \lim_{\rho\to 0} \frac{ D\car{E}(B_\rho(0))}{|D\car{E}|(B_\rho(0))|}
\end{equation}
and, without loss of generality, we assume that it coincides with
the vector $e_d$ corresponding to the last coordinate of $y\in \mathbb{R}^d$.

Also, if we let $E_\rho:= E/\rho$, the sets $E_\rho$, $E_\rho^c$,
$\partial E_\rho$
converge in $B_1(0)$, in the Hausdorff sense (thanks to the uniform
density estimates), respectively to $\{y_d\ge 0\}$, $\{y_d=0\}$, $\{y_d\le 0\}$.
We also let $z_\rho(y) := z(\rho y)$ and $g_\rho(y) := g(\rho y)$,
in particular $-\textup{div}\, z_\rho=\rho g_\rho$. We let
 \begin{equation}\label{defomega}
 \omega(\rho)\ =\  \sup_{x\in\Omega}\|g\|_{L^d(B_\rho(x)\cap\Omega)}\,
 \end{equation}
which is continuously increasing and goes to $0$ as $\rho\to 0$,
since $|g|^d$ is equi-integrable.

We introduce the following notation: a point in $\mathbb{R}^d$ is denoted by
$y=(y',y_d)$, with $y'\in \mathbb{R}^{d-1}$. We let $D_s:=\{|y'|\le s\}$, $\bar z:= \nabla F(\ov{\nu})$ and
$D_s^t = \{ D_s+ \lambda \bar z: |\lambda|\le t\}$ and denote with $\partial D_s$ the relative boundary
of $D_s$ in $\{y_d=0\}$.

We choose $s\le 1$, $0< t\le s$, ($t$ is chosen small enough
so that $D_s^t\subset B_1(0)$, that is $t < (1/|\bar z|)\sqrt{1-s^2}$).
We integrate in $D_s^t$ the divergence $\rho g_\rho=-\textup{div}\, z_\rho=\textup{div}\,(\bar z-z_\rho)$
against the function $(2\car{E}-1)t-\frac{ \ov{\nu}\cdot y}{F(\ov{\nu})} $, which vanishes
for $y_d=\pm t F(\ov{\nu}) $ if $\rho$ is small enough (given $t>0$), so
that $\partial E_\rho\cap B_1(0)\subset \{|y_d|\le t F(\ov{\nu})\}$. For $y$ on the lateral boundary of the cylinder $D_s^t$, let $\xi(y)$ be the internal normal to $\partial D_s+(-t,t)\bar z$ at the point $y$.
Using the fact that $z_\rho$ is a calibration for $E_\rho$,
we easily get that for almost all $s$,
\begin{multline}\label{bigone}
\int_{D_s^t} \rho g_\rho\left((2\car{E}-1)t-\frac{ \ov{\nu}\cdot y}{F(\ov{\nu})}\right)\,dy
\\ =\
\int_{\partial D_s+(-t,t)\bar z} \left((2\car{E}-1)t-\frac{ \ov{\nu}\cdot y}{F(\ov{\nu})}\right)
[(\bar z-z_\rho),\xi]\,d\mathcal{H}^{d-1}
\\-\,2t\int_{\partial E_\rho\cap D_s^t}
\left(\bar z\cdot \nu^{E_\rho}-F(\nu^{E_\rho})\right)\,d\mathcal{H}^{d-1}
\,+\, \int_{D_s^t} \left(1-\frac{z_\rho\cdot\ov{\nu}}{F(\ov{\nu})}\right)\,dy\,.
\end{multline}
Now since $F^\circ(\nabla F(\ov{\nu}))=1$, there holds $\bar z\cdot \nu^{E_\rho}-F(\nu^{E_\rho})\le 0$ and using that $\bar z\cdot \xi(y)=0$ on $\partial D_s+(-t,t)\bar z$, we get
\begin{multline}\label{bigtwo}
\int_{D_s^t} \left(1-\frac{z_\rho\cdot\ov{\nu}}{F(\bar \nu)}\right)\,dy\le\
\int_{D_s^t} \rho g_\rho\left((2\car{E}-1)t-\frac{ \ov{\nu}\cdot y}{F(\ov{\nu})}\right)\,dy
\\\int_{\partial D_s+(-t,t)\bar z} \left((2\car{E}-1)t-\frac{ \ov{\nu}\cdot y}{F(\ov{\nu})}\right)
z_\rho\cdot \xi\,d\mathcal{H}^{d-1}\,.
\end{multline}

We claim that for $|\xi|\le 1$ with $\xi\cdot \bar z=0$, there holds
\begin{equation}\label{star}
 (\xi\cdot z_\rho)^2\le C(F(\ov{\nu})-\ov{\nu}\cdot z_\rho)
\end{equation}
 Since
\[(\xi\cdot z_\rho)^2\le |z_\rho|^2-[z_\rho\cdot (\bar z/|\bar z|)]^2\]
it is enough to prove
\[|z_\rho|^2-[z_\rho\cdot (\bar z/|\bar z|)]^2\le  C(F(\ov{\nu})-\ov{\nu}\cdot z_\rho).\]

Using that $\ov{\nu}/F(\ov{\nu})=\nabla F^\circ (\bar z)$, from \eqref{strictconvF} applied to $F^{\circ}$ together with $F^\circ(\bar z)=1\ge F^{\circ}(z_\rho)$, we find
\[(F(\ov{\nu})-\ov{\nu}\cdot z_\rho)=F(\ov{\nu}) (1-z_\rho\cdot \nabla F^\circ(\bar z))\ge C|z_\rho-\bar z|^2.\]
which readily implies  \eqref{star}. We thus have
\begin{multline}\label{big4}
\int_{\partial D_s+(-t,t)\bar z}
\left((2\car{E_\rho} -1) t -\frac{\ov{\nu}\cdot y}{F(\ov{\nu})} \right)
\,\left(z_\rho\cdot \xi\right)\,d\mathcal{H}^{d-1}
\\ \le\,
2 C\sqrt{F(\ov{\nu})} t \int_{\partial D_s+(-t,t)\bar z}
\sqrt{1-\frac{z_\rho\cdot\ov{\nu}}{F(\ov{\nu})}}\;d\mathcal{H}^{d-1}
\\ \le\,
2CF(\ov{\nu}) t\sqrt{t} \left(
\int_{\partial D_s+(-t,t)\bar z}
 \left(1-\frac{z_\rho\cdot\ov{\nu}}{F(\ov{\nu})}\right) \,d\mathcal{H}^{d-1}
\right)^\frac{1}{2} \sqrt{\mathcal{H}^{d-2}(\partial D_s)}\,.
\end{multline}
Now, we also have
\begin{multline}\label{big5}
\rho\int_{D_s^t}\left((2\car{E_\rho} -1) t -\frac{\ov{\nu}\cdot y}{F(\ov{\nu})}\right)
g_\rho\,dy\,\le\,2t\rho^{1-d}\int_{D_{\rho s}^{\rho t}}g \,dy
\\
\le\, 2t \rho^{1-d}\|g\|_{L^d(B_{\rho s}(0))}  |D_{\rho s}^{\rho t}|^{1-1/d}
\,\le\, c t^{2-1/d} s^{d-2+1/d} \omega(\rho s)
\end{multline}
where here, $c=2\mathcal{H}^{d-1}(D_1)^{1-1/d}$,
and $\omega$ is defined in~\eqref{defomega}.

We choose $a<1$, close to $1$, and  $t\in (0,(1/|\bar z|)\sqrt{1-a^2})$.
If $\rho>0$ is small enough
(so that $\partial E_\rho\cap B_1$ is in $\{ |y_d|\le t F(\ov{\nu})\}$),
letting $f(s):=\int_{D^t_s} \left(1-\frac{z_\rho\cdot\ov{\nu}}{F(\ov{\nu})}\right)dy$,
we deduce from~\eqref{bigtwo}, \eqref{big4} and~\eqref{big5}
that for a.e. $s$ with $t\le s\le a$, one has (possibly increasing
the constant $c$)
\begin{equation}\label{finalestimate}
f(s)^2\ \le\  c\left( s^{d-2} t^3  f'(s)
\,+\,t^{4-2/d} s^{2d-4+2/d} \omega(\rho s)^2\right)\,.
\end{equation}
Unfortunately, this estimate does not seem to  give much information for $d>3$.
It seems it allows to conclude only whenever $d\in \{2,3\}$.
Since the case $d=2$
is simpler, we focus on $d=3$. Estimate~\eqref{finalestimate} becomes
\begin{equation}\label{finalestimate3D}
f(s)^2\ \le\  c\left( s t^3  f'(s)
\,+\,  t^{10/3} s^{8/3} \omega(\rho s)^2\right)\,.
\end{equation}
Given $M>0$, we fix a value $t>0$ such that $\log(a/t)\ge c M$.
If $\rho$ is chosen small enough, then $\partial E_\rho\cap B_1(0)
\subset \{ |y_d| < t F(\ov{\nu})\}$, and~\eqref{finalestimate3D} holds.
It yields (assuming $f(t)>0$, but if not, then the proposition
is proved)
\begin{equation}\label{ff}
-\frac{f'(s)}{f(s)^2}\,+\,\frac{1}{c t^3} \frac{1}{s} \ \le\
ct^{1/3} s^{5/3} \frac{\omega(\rho s)^2}{f(s)^2}
\ \le\
ct^{1/3} s^{5/3} \frac{\omega(a \rho )^2}{f(t)^2}
\end{equation}
where we have used the fact that $t\le s\le a$ and  $f,\omega$
are nondecreasing.
Integrating \eqref{ff} from $t$ to $a$, after multiplication by $t^3$ we obtain
\[
\frac{t^3}{f(a)} - \frac{t^3}{f(t)} + \frac{\log(a/t)}{c}
\ \le\ \frac{3c}{8} t^{10/3}(a^{8/3}-t^{8/3})\frac{\omega(a\rho)^2}{f(t)^2}\,.
\]
Hence we get
\begin{equation}\label{finalestimatef}
\frac{t^3}{f(t)} + c a^{8/3} t^{-8/3}\omega(a\rho)^2 \frac{t^6}{f(t)^2} \ \ge\ M.
\end{equation}

Eventually, we observe that
\[
f(t)=\int_{D_t^t}\left(1-\frac{z(\rho y)\cdot\ov{\nu}}{F(\ov{\nu})}\right)\,dy = \frac{1}{\rho^d}\int_{D_{\rho t}^{\rho t}} \left(1-\frac{z(x)\cdot\ov{\nu}}{F(\ov{\nu})}\right)\,dx\,,
\]
so that \eqref{finalestimatef} can be rewritten
\begin{equation}\label{almostdone}
\left(\frac{\displaystyle \int_{D_{\rho t}^{\rho t}}
\left(1-\frac{z(x)\cdot\ov{\nu}}{F(\ov{\nu})}\right)dx}{(\rho t)^3}\right)^{-1}
\ \ge\
\frac{-1+\sqrt{ 1+4 M c a^{8/3} t^{-8/3} \omega(a\rho)^2}}
{2c a^{8/3}t^{-8/3} \omega(a\rho)^2}
\end{equation}
The value of $t$ being fixed, we can choose the value of $\rho$ small enough
in order to have
$4 M c a^{8/3} t^{-8/3} \omega(a\rho)^2 < 1$, and (using
$\sqrt{1+X}\ge 1+X/2 - X^2/8$ if $X\in (0,1)$), \eqref{almostdone}
yields
\begin{equation}\label{almostdone2}
\left(\frac{\displaystyle \int_{D_{\rho t}^{\rho t}}
\left(1-\frac{z\cdot\ov{\nu}}{F(\ov{\nu})}\right)\, dy}{(\rho t)^3}\right)^{-1}
\ \ge\  M - M^2 ca^{8/3} t^{-8/3}\omega(a\rho)^2\ \ge\ \frac{3}{4}M\,.
\end{equation}
It  follows that
\begin{equation}
\limsup_{\varepsilon\to 0} \frac{\displaystyle \int_{D_{\varepsilon}^{\varepsilon}} \left(1-\frac{z\cdot\ov{\nu}}{F(\ov{\nu})}\right)\,dy}{\varepsilon^3}
\ \le\ \frac{4}{3}M^{-1}
\end{equation}
and since $M$ is arbitrary, $0$ is indeed a Lebesgue point of $z$, with
value $\bar z=\nabla F(\ov{\nu})$ (recall that $1-\frac{z(x)\cdot\ov{\nu}}{F(\ov{\nu})} \ge (C/F(\ov{\nu})) |z(x)-\bar z|^2$).

{\it Step 2.} When $F$ depends also on the $x$ variable, the proof follows along the same lines
as in {\it Step 1}, taking into account the errors terms in \eqref{bigtwo} and \eqref{big4}.
Keeping the same notations as in {\it Step 1}
and setting $\bar z:= \nabla_p F(0,\bar \nu)$ we find that since  $F^\circ(0,\bar z)\le 1$, there holds $\bar z\cdot \nu^{E_\rho}\le F(0,\nu^{E_\rho})$ and thus
\[\int_{\partial E_\rho \cap D_s^t} \bar z\cdot \nu^{E_\rho}-F(\rho x, \nu^{E_\rho}) d\mathcal{H}^{d-1}\le\int_{\partial E_\rho \cap D_s^t} |F(0,\nu^{E_\rho})-F(\rho x,\nu^{E_\rho})| d\mathcal{H}^{d-1}
\le C \rho s^{d-1}
\]
where the last inequality follows from $t\le s$ and the minimality of $E_\rho$ inside $D_s^t$. Now since
\[\left( F^\circ\right)^2(0,z_\rho)-\left( F^\circ\right)^2(\rho x,z_\rho)\ge \left( F^\circ\right)^2(0,z_\rho)-1\ge 2 \frac{\bar \nu}{F(0,\bar \nu)} \cdot (z_\rho-z) +\delta^2|z_\rho-z|^2\]
we find that \eqref{star} transforms into,
\[(\xi\cdot z_\rho)^2\le C \left[(F(0,\bar \nu)-\bar \nu\cdot z_\rho) +(\left( F^\circ\right)^2(0,z_\rho)-\left( F^\circ\right)^2(\rho x,z_\rho))\right]\]
for every $|\xi|\le 1$ and $\xi\cdot \bar z=0$, from which we get
\begin{multline*}
 \int_{\partial D_s+(-t,t)\bar z}
\left((2\car{E_\rho} -1) t -\frac{\ov{\nu}\cdot y}{F(\ov{\nu})} \right)
\,\left(z_\rho\cdot \xi\right)\,d\mathcal{H}^{d-1}
\\ \le\,
2CF(0,\ov{\nu}) t\sqrt{t} \left(
\int_{\partial D_s+(-t,t)\bar z}
 \left(1-\frac{z_\rho\cdot\ov{\nu}}{F(0,\ov{\nu})}\right) \,d\mathcal{H}^{d-1}
\right)^\frac{1}{2} \sqrt{\mathcal{H}^{d-2}(\partial D_s)}\\
+ 2Ct \,\int_{\partial D_s+(-t,t)\bar z}
 \left|\left( F^\circ\right)^2(0,z_\rho)-\left( F^\circ\right)^2(\rho x,z_\rho)\right|^{1/2}\, d\mathcal{H}^{d-1}\\
\le \,C F(0,\ov{\nu}) t\sqrt{t} \left(
\int_{\partial D_s+(-t,t)\bar z}
 \left(1-\frac{z_\rho\cdot\ov{\nu}}{F(0,\ov{\nu})}\right) \,d\mathcal{H}^{d-1}
\right)^\frac{1}{2} \sqrt{\mathcal{H}^{d-2}(\partial D_s)}+ Ct \rho^{1/2} s^{d-1} t\,.
\end{multline*}
Using these estimates, we finally get that, setting as before $f(s):=\int_{D^t_s} \left(1-\frac{z_\rho\cdot\ov{\nu}}{F(0,\ov{\nu})}\right)dy$, there holds
\[f(s)^2\ \le\  c\left( s^{d-2} t^3  f'(s)
\,+\,t^{4-2/d} s^{2d-4+2/d} \omega(\rho s)^2+ \rho t s^{d-1}+\rho^{1/2} t^2 s^{d-1}\right)\,.\]
{}From this inequality, the proof can be concluded exactly as in {\it Step 1}.
\end{proof}

\begin{remark}\label{remDu}\rm
Assuming $F$ has some regularity (Lipschitz in the first variable,
and $C^{2,\beta}$ and even in the second, see~\cite{SchoenSimon}), then for
$d=2$ or $d=3$ and $g\in L^p(\Omega)$ with $p>d$, $\partial E$ is of class $C^{1,\alpha}$ for some $\alpha>0$. In this case, \eqref{euleru} holds everywhere in $\textup{spt}\,(Du)$.
\end{remark}

Eventually, we can also give a locally uniform convergence result (valid
in dimension $d=2,3$, with the assumption\footnote{Probably just technical.} that $F$ is even in dimension 3).
\begin{proposition}
For all $x\in\Omega$ we let $$z_\rho(x):=\frac{1}{|B_\rho(0)|}\int_{B_\rho(x)\cap\Omega}z \,dy\,.$$
Then, $F^\circ(x,z_\rho(x))\to 1$ locally uniformly on $S$.
\end{proposition}
\begin{proof}
Given $K\subset \Omega$ a compact set, we can check that
for any $t>0$, there exists $\rho_0>0$ such that for any $x\in K\cap S$,
if $E^x$ is the level set of $u$ through $x$, then
for any $\rho\le\rho_0$, the boundary of
$(E^x-x)/\rho\cap B_1(0)$ lies in a strip of width $2t$, that is, there is
$\ov{\nu}^x_\rho\in \mathbb{S}^{d-1}$ with $\partial((E^x-x)/\rho)\cap B_1(0)\subset
\{|y\cdot\ov{\nu}^x_\rho |\le t\}$).

Indeed, if this is not the case, one can find $t>0$, $\rho_k\to 0$, $x_k\in K\cap S$,
such that $\partial((E^{x_k}-x_k)/\rho_k)\cap B_1(0)$ is not contained
in any strip of width $2t$. Up to a subsequence we may assume
that $x_k\to x\in K\cap S$, and from the bound on the perimeter, that $(E^{x_k}-x_k)/\rho_k$ converges (in the Kuratowski sense) to a local minimizer of $\int_{\partial E} F(0,\nu^E) d\mathcal{H}^{d-1}$
 and is thus a halfspace.\footnote{If $d=2$, this Bernstein result readily follows from the strict convexity of $F$, see \cite[Prop 3.6]{CGN}
whereas for $d=3$, see \cite[Thm.~4.1]{white}, where it is assumed that $F$ is even.
In the case of the area i.e when $F(x,Du)=|Du|$ and $d\le 7$, see also \cite[Rem 3.2]{GoMa}.} It yields that $\partial((E^{x_k}-x_k)/\rho_k)\cap B_1(0)$
converges in the Hausdorff sense (thanks to the density estimates) to a
hyperplane. We easily obtain a contradiction.

The thesis follows when we observe that the proof of Proposition~\ref{propLebesgue}
can be reproduced by replacing
the direction $\nu^{E^x}(x)$ (which exists only if $x$ lies in the reduced boundary of $E^x$)
with the direction $\ov{\nu}^x_\rho$ given above.
\end{proof}
\begin{remark}\rm
In the Euclidean case ($F=|\,.\,|$) it has already been observed in~\cite{AmbPao}
that the blow-ups are flat at each point of the boundary of a set
with curvature in $L^d$ (but for a closed set of maximal dimension $d-8$),
however the spiral example in~\cite{GonzalezMassariTamanini} shows that
even if $d=2$, the orientation of the limit line may not be unique.
\end{remark}
\subsection{A counterexample.}
We provide an example where $g\in L^{d-\varepsilon}(\Omega)$, with $\varepsilon>0$ arbitrarily small,
and Theorem \ref{propLebesgue} does not hold.

Let $\Omega=B_1(0)$ be the unit ball of $\mathbb{R}^d$ and let $E=\Omega\cap\{ x_d\le 0\}$. We shall construct a vector
field $z:\Omega\to \mathbb{R}^d$ such that $z=\nu^E$ on $\partial E\cap \Omega$, $|z|\le 1$ everywhere in $\Omega$, ${\rm div}z\in L^{d-\varepsilon}(\Omega)$,
but $0$ is not a Lebesgue point of $z$. Notice that $E$ minimizes the functional \eqref{geomper} with $g={\rm div}z$.

Letting $r_n\to 0$ be a decreasing sequence to be determined later, and let $B_n=B_{r_n}(x_n)$ with $x_n=2 r_n e_d$.
Without loss of generality, we may assume $r_{n+1}<r_n/4$ so that the balls $B_n$ are all disjoint.
We define the vector field $z$ as follows: $z(x)=e_d$ if $x\in\Omega\setminus\cup_n B_n$, and
$z(x)=|x-x_n|e_d$ if $x\in B_n$.
It follows that ${\rm div}z=0$ in $\Omega\setminus\cup_n B_n$ and $|{\rm div}z|\le 1/r_n$ in $B_n$, so that
\[
\int_\Omega |{\rm div}z|^{d-\varepsilon} \,dx = \sum_n \int_{B_n} |{\rm div}z|^{d-\varepsilon} \,dx \le \omega_d \sum_n r_n^\varepsilon<+\infty
\]
if we choose $r_n$ converging to zero sufficiently fast, so that $g=-\textup{div}\, z\in L^{d-\varepsilon}(\Omega)$.

However, since $z\cdot e_d\le 1/2$ in $B_{r_n/2}(x_n)$, we also have
\[
\int_{B_{3r_n}(0)}z\cdot e_d\,dx \le \left|B_{3r_n}(0)\right| -\frac 12 \left|B_{r_n/2}(x_n)\right|
\]
so that
\[
\frac{1}{\left|B_{3r_n}(0)\right|}\, \int_{B_{3r_n}(0)}z\cdot e_d\,dx
\le 1-\frac{1}{6^d}<1\,.
\]
On the other hand, for $\delta\in (0,1/6^d)$ we have
\[
\frac{1}{\left|B_{r_n}(0)\right|}\,
\int_{B_{r_n}(0)}z\cdot e_d\,dx \ge
\frac{1}{\left|B_{r_n}(0)\right|}\left(
\left|B_{r_n}(0)\right| - \sum_{i=n+1}^\infty \left|B_{r_i}(x_i)\right|\right)\ge 1-\delta\,,
\]
if we take the sequence $r_n$ converging to $0$ sufficiently fast.
It follows that $0$ is not a Lebesgue point of $z$.