# Conformal deformations of initial data sets to the strict dominant energy
condition using a spacetime Poisson equation

Jaroslaw S. Jaracz

###### Abstract.

We give an alternate proof of one of the results given in [Lee2022] showing
that initial data sets with boundary for the Einstein equations $(M,g,k)$
satisfying the dominant energy condition can be conformally deformed to the
strict dominant energy condition, while preserving the character of the
boundary (minimal, future trapped, or past trapped) while changing the area of
the boundary and ADM energy of the initial data set by an arbitrarily small
amount. The proof relies on solving an equation that looks like the equation
for spacetime harmonic functions studied in [HirschKazarasKhuri], but with a
Neumann boundary condition and non-zero right hand side, which we refer to as
a spacetime Poisson equation. One advantage of this method of proof is that
the conformal deformation is explicitly constructed as a solution to a PDE, as
opposed to only knowing the solution exists via an application of the implicit
function theorem as in [Lee2022]. We restrict ourselves to the physically
relevant case of a $3$-manifold $M$, though the proof can be generalized to
higher dimensions.

##  1\. Introduction

When studying initial data sets $(M,g,k)$ for the Einstein equations, one is
usually interested in initial data sets satisfying certain energy conditions,
such as the dominant energy condition

$\displaystyle\mu\geq|J|_{g}$

(see Section 2.2 for definitions of the given quantities) which can be
interpreted as saying that the speed of energy flow never exceeds that of the
speed of light. However, for various technical reasons, it is easier to deal
with the strict dominant energy condition

$\displaystyle\mu>|J|_{g}.$

The strict condition gives enough wiggle room to make certain technical
arguments a lot easier. Therefore, when proving theorems for which the
dominant energy condition is valid, one would like to to say it is possible to
approximate the initial data sets by ones satisfying the strict condition in
some appropriate sense.

The first example of this sort of argument was given in [SchoenYau2], and the
deformation to the strict energy condition was claimed in Lemma 1 of that
paper, using the implicit function theorem. However, as was pointed out in
[EichmairHuangLeeSchoen2016], there was actually a flaw in this argument as
the implicit function theorem was not applicable in that case.

It was only in [EichmairHuangLeeSchoen2016] that a correct argument was given
(for manifolds without boundary), by introducing something called the modified
Lie derivative, which allowed for the correct application of the implicit
function theorem on Banach spaces. This was then generalized to the case of
manifolds with boundary in [Lee2022]. Now there are two things regarding these
proofs which are not ideal. One is that the introduction of the modified Lie
derivative and subsequent calculations can be somewhat tricky to follow along.
The second is that since the proofs rely on the implicit function theorem,
they only say that a deformation exists, but not how to obtain it.

While studying conformal deformations for applications to the Penrose
inequality (see Section 4), we stumbled across an alternative proof of such
deformations, which does not rely on the implicit function theorem. On a
conceptual level, the deformation is given as a solution to a certain quasi-
linear elliptic PDE with Neumann boundary conditions. This also potentially
has some useful practical applications. For example, if one wants to
numerically study some initial data set, and wants to make sure the initial
data set satisfies the strict dominant energy condition, one can explicitly
calculate the solution to the aforementioned PDE. We thought this was
sufficiently interesting from both the theoretical and practical point of view
to be worth writing up.

We do mention one advantage of the proofs in [EichmairHuangLeeSchoen2016,
Lee2022] is that they yield so-called harmonic asymptotics (which are in
themselves useful), and it is unclear if our proof can be modified to yield
these asymptotics as well. They also give results for initial data sets with
lower regularity, though our proof should generalize to this case in a
straightforward way.

We focus on the physically relevant case where $M$ is $3$-dimensional, though
the ideas can be generalized to arbitrary dimension. To simplify presentation,
we will assume that $M$ has asymptotially flat coordinates (see Section 2.1)
in which

$\displaystyle\begin{split}M=\mathbb{R}^{3}\setminus U\\\ 0\notin
M\end{split}$ (1.1)

where $U$ is some domain in $\mathbb{R}^{3}$ with possibly several components,
with $\partial U=\partial M$ smooth and compact with finitely many components.
By choosing the coordinates so that $0\notin M$ we have $r(x)\geq c_{0}>0$ for
all $x\in M$ and so we can take some constant $C$ such that (2.1) holds for
all $x\in M$.

This is done for convenience and to make the idea of the proof more
transparent. The case of primary interest is that of one asymptotically flat
end where $\partial M$ is nonempty and an outermost apparent horizon (defined
below). In that case using a change of coordinates we can always choose
coordinates such that (1.1) holds.

The theorem still applies to $M$ not of the form (1.1) but in that case the
function $r(x)$ appearing in (1.2) has to be replaced by some appropriately
defined weight function $\sigma(x)$. This $\sigma(x)$ could be taken to be
some smooth function defined to satisfy $\sigma\equiv 1$ on some compact $K$
and $\sigma(x)=r(x)$ sufficiently far in the asymptotically flat end. All the
subsequent proofs can be generalized to this case in a straightforward, though
slightly messy manner.

The main theorem we prove in this paper is:

######  Theorem 1.1.

Let $(M,g,k)$ be a smooth asymptotically flat initial data set for the
Einstein equations with $M$ of the form (1.1), satisfying the dominant energy
condition, with one asymptotically flat end and compact boundary $\partial M$
of area $A$, satisfying the fall-off conditions (2.1) with constant $C$, and
with ADM energy $E$. Then there exists a $\psi\in C^{2,\alpha}_{-1}({M})$
solving the Neumann problem

$\displaystyle\begin{split}\Delta_{g}\psi+\frac{C}{r^{2}}|\nabla\psi|_{g}&=-\frac{1}{r^{4}}\\\
\left.\frac{\partial\psi}{\partial\nu}\right|_{\partial M}&=0.\end{split}$
(1.2)

Moreover, if we define $u_{\epsilon}=1+\epsilon\psi$ for all sufficiently
small $\epsilon>0$ then the initial data set
$(M,g_{\epsilon},k_{\epsilon})=(M,u_{\epsilon}^{4}g,u_{\epsilon}^{2}k)$ is
asymptotically flat, satisfies the strict dominant energy condition, and there
exists some constant $\mathcal{C}>0$ such that

$\displaystyle|E-E_{\epsilon}|\leq\mathcal{C}\epsilon,\quad|A-A_{\epsilon}|\leq\mathcal{C}\epsilon$
(1.3)

where $E_{\epsilon}$ and $A_{\epsilon}$ are the ADM energy and boundary area
of $g_{\epsilon}$. Moreover, the sign of the null expansions at the boundary
is the same for $g$ and $g_{\epsilon}$, so in particular apparent horizon
boundaries remain apparent horizons.

We mention that the right hand side of (1.2) can be replaced by any negative
function which vanishes sufficiently fast at infinity, giving further
(potentially useful) control over the deformation. The proof closely follows
[HirschKazarasKhuri] where a similar equation, but with Dirichlet boundary
condition, is studied. The proof consists of the application of standard
tools, such as the Leray-Schauder fixed point theorem, maximum principles, and
some simple elliptic estimates.

Finally, we mention that a linear elliptic equation with Neumann boundary
condition was used to study the scalar curvature of initial data sets $(M,g)$
in [Lee2023].

##  2\. Background Material

###  2.1. Asymptotic Flatness and ADM Energy

We consider an asymptotically flat initial data set $(M,g,k)$ where $M$ is a
$3$-manifold, $g$ a Riemannian metric, and $k$ is a symmetric $2$-tensor, the
extrinsic curvature. We assume the manifold is smooth and has a single
asymptotically flat end and some smooth compact boundary $\partial M$ with
possibly multiple components. These components are not necessarily apparent
horizons, though the apparent horizon case is the one of greatest interest. We
assume that $g$ and $k$ are of class $C^{\infty}({M})$.

An initial data set having a single asymptotically flat end means that there
is some compact set $K$ such that $M\setminus K$ is diffeomorphic to the
complement of a ball in $\mathbb{R}^{3}$. This diffeomorphic region is
referred to as the asymptotically flat end. It is assumed that on this
asymptotically flat end $g$ and $k$ satisfy certain fall-off conditions. In
our case, we take these to be

$\displaystyle|D^{\lambda}(g_{ij}-\delta_{ij})|\leq
Cr^{-1-|\lambda|},\quad|R|\leq Cr^{-3},\quad|D^{\lambda}k|\leq
Cr^{-2-|\lambda|}$ (2.1)

for $|\lambda|\leq 2$. Here, $C$ is some constant, $R$ is the scalar
curvature, $\delta$ is the Euclidean metric, $r=\sqrt{x^{2}+y^{2}+z^{2}}$ the
standard Euclidean radius, $D^{\lambda}$ is a derivative operator with respect
to the Euclidean coordinates, and $\lambda$ is a multi-index. We have

$|k|^{2}=k_{ij}k^{ij},\quad Tr_{g}k=g^{ij}k_{ij}$

as usual, where the indices are raised with $g$.

For an asymptotically flat end, the ADM energy is defined by

$\displaystyle
E_{ADM}=\lim_{r\rightarrow\infty}\frac{1}{16\pi}\sum_{i,j}\int_{S_{r}}(g_{ij,i}-g_{ii,j})\nu^{j}dS_{r}$
(2.2)

where $S_{r}$ are coordinate spheres of radius $r$ and $\nu^{j}$ is the
outward unit normal [ADM].

###  2.2. Einstein Constraint Equations and Energy Conditions

The Einstein constraint equations for an initial data set $(M,g,k)$ are

$\displaystyle\begin{split}\mu&=\frac{1}{2}\left(R+(Tr_{g}k)^{2}-|k|_{g}^{2}\right)\\\
J_{i}&=\nabla^{j}(k_{ij}-(Tr_{g}k)g_{ij})\end{split}$ (2.3)

where $R$ is the scalar curvature, $\nabla^{j}$ denotes covariant
differentiation, $\mu$ is the mass-energy density, and $J_{i}$ the components
of the momentum density. Then the dominant energy condition (DEC) is that the
inequality

$\mu\geq|J|_{g}$ (2.4)

is satisfied everywhere in the manifold.

###  2.3. Null expansions and apparent horizons

Given a (two sided) codimension-$1$ hypersurface in $M$ one can calculate the
future $(+)$ and past $(-)$ null expansions at each point of the surface,
which are defined by

$\displaystyle\theta_{\pm}=H\pm Tr_{S}k.$ (2.5)

Here $H$ is the mean curvature computed with respect to the unit normal $\nu$
pointing towards spatial infinity in a chosen AF end and $Tr_{S}k$ denotes the
trace of $k$ restricted to the surface $S$ computed with respect to the
induced metric on the surface.

The futre and past null expansions measure the convergence and divergence of
past and future directed null geodesics. Surfaces with $\theta_{\pm}<0$ are
referred to as future and past trapped. Surfaces with $\theta_{\pm}$ are
referred to as future and past apparent horizons and are a popular way of
modelling the surface of a black hole in an initial data set.

###  2.4. Conformal Change

A conformal change of an initial data set $(M,g,k)$ to a new initial data set
$(M,\bar{g},\bar{k})$ is defined by

$\displaystyle(M,g,k)\rightarrow(M,\overline{g},\overline{k})=(M,u^{4}g,u^{2}k).$

Here $u>0$ is some non-negative function on $M$. It is then useful to know how
the relevant quantities transform under this conformal change. One has

$\displaystyle\begin{split}\overline{\mu}&=u^{-4}\left(\mu-4u^{-1}\Delta_{g}u\right)\\\
\overline{J}_{i}&=u^{-2}\left(J_{i}+4u^{-1}k_{ij}\nabla^{j}u\right).\end{split}$
(2.6)

See the formulas in Lemma 1 of [SchoenYau2] where our $\overline{J}_{i}$ is
related to the $K_{i}$ defined there by $\overline{J}_{i}=u^{-2}K_{i}$. Though
the proof of the lemma incorrectly applies the implicit function theorem on
Banach spaces, the formulas for conformal change are correct.

We want to derive some condition on $u$ such that the DEC is preserved. Notice
that for a $1$-form $X_{i}$ we have

$\displaystyle|X_{i}|_{\overline{g}}^{2}=\overline{g}^{ij}X_{i}X_{j}=u^{-4}g^{ij}X_{i}X_{j}=u^{-4}|X_{i}|_{g}^{2}$
(2.7)

so that

$\displaystyle|X_{i}|_{\overline{g}}=u^{-2}|X_{i}|_{g}.$ (2.8)

Combining this with (2.6) and the triangle inequality we obtain

$\displaystyle|\overline{J}_{i}|_{\overline{g}}\leq|u^{-2}J_{i}|_{\overline{g}}+|4u^{-3}k_{ij}\nabla^{j}u|_{\overline{g}}=u^{-4}|J_{i}|_{{g}}+4u^{-5}|k_{ij}\nabla^{j}u|_{{g}}$
(2.9)

Therefore, if we have

$\displaystyle\Delta_{g}u\leq-|k_{ij}\nabla^{j}u|_{g}$ (2.10)

then

$\displaystyle\overline{\mu}$ $\displaystyle=u^{-4}\mu-4u^{-5}\Delta_{g}u\geq
u^{-4}\mu+4u^{-5}|k_{ij}\nabla^{j}u|_{g}\geq
u^{-4}|J_{i}|_{{g}}+4u^{-5}|k_{ij}\nabla^{j}u|_{{g}}\geq|\overline{J}_{i}|_{\overline{g}}$
(2.11)

and the dominant energy condition is preserved. Of course, then the most
natural equation to try to solve is

$\displaystyle\Delta_{g}u=-|k_{ij}\nabla^{j}u|_{g}.$ (2.12)

However, since using (2.1) we can estimate

$\displaystyle|k_{ij}\nabla^{j}u|_{g}\leq\frac{C}{r^{2}}|\nabla u|_{g}$ (2.13)

then if we look to solve

$\displaystyle\Delta_{g}u=-\frac{C}{r^{2}}|\nabla u|_{g}$ (2.14)

then we get

$\displaystyle\Delta_{g}u=-\frac{C}{r^{2}}|\nabla
u|_{g}\leq-|k_{ij}\nabla^{j}u|_{g}$

preserving the dominant energy condition. Notice that (2.14) is of the same
form as the equation for spacetime harmonic functions for which existence
results were established in [HirschKazarasKhuri], which we will exploit. The
equation studied in [HirschKazarasKhuri] was of the form

$\Delta_{g}u+G(x)|\nabla u|_{g}=0$ (2.15)

and since the solutions were referred to as spacetime harmonic functions, it
makes sense to refer to (2.15) as a spacetime Laplace equation. Therefore it
also makes sense to refer to an equation of the form

$\displaystyle\Delta_{g}u+G(x)|\nabla u|_{g}=F(x)$ (2.16)

as a spacetime Poisson equation, hence our terminology.

Let $S$ be a two sided codimension-$1$ surface in $M$ with mean curvature $H$
computed with respect to $g$ with the chosen unit normal $\nu$. Then the mean
curvature computed with respect to $\bar{g}$ with respect to the unit normal
$\overline{\nu}=u^{-2}\nu$ is given by the standard formula

$\displaystyle\overline{H}=u^{-2}\left(H+4u^{-1}\langle\nabla
u,\nu\rangle_{g}\right),$ (2.17)

see formula (7) in [EichmairHuangLeeSchoen2016].

If we denote by $h$ the induced metric on $S$ then $\overline{h}=u^{4}h$. If
we let $K=k_{S}$ denote $k$ restricted to $S$, and similarly
$\overline{K}=\overline{k}_{S}$ then we get

$\displaystyle\overline{K}=u^{2}K$ (2.18)

so that

$\displaystyle
Tr_{S}\overline{k}=\overline{h}^{ij}\overline{K}_{ij}=u^{-4}h^{ij}u^{2}K_{ij}=u^{-2}h^{ij}K_{ij}=u^{-2}Tr_{S}k.$
(2.19)

Therefore we have

$\displaystyle\overline{\theta}_{\pm}=u^{-2}\theta_{\pm}+4u^{-3}\langle\nabla
u,\nu\rangle_{g}$ (2.20)

by combining the above formulas.

##  3\. Proof of Theorem 1.1

We first prove all the claims of the theorem assuming a solution to (1.2)
exists. We then prove the existence of $\psi$. For the definitions of the
weighted Hölder spaces $C^{k,\alpha}_{-\beta}$ we refer the reader to
[EichmairHuangLeeSchoen2016, Chaljub1979, Weinstein_2005]. Suppose a solution
$\psi\in C^{2,\alpha}_{-1}({M})$ to the problem (1.2) exists. Then in
particular $\psi$ is bounded. Define

$\displaystyle u_{\epsilon}=1+\epsilon\psi.$ (3.1)

so for all sufficiently small $\epsilon>0$ we have that $1+\epsilon\psi>0$ so
$g_{\epsilon}=u_{\epsilon}^{4}g$ is a Riemannian metric. Due to the fact that
$\psi\in C^{2,\alpha}_{-1}({M})$ the function $u_{\epsilon}$ tends towards $1$
as $|x|\rightarrow\infty$ and the derivatives $\partial_{i}u_{\epsilon}$ and
$\partial_{ij}u_{\epsilon}$ vanish as $O(r^{-2})$ and $O(r^{-3})$
respectively. Therefore, also defining $k_{\epsilon}=u_{\epsilon}^{2}k$ yields
that $(M,g_{\epsilon},k_{\epsilon})$ is an asymptotically flat initial data
set.

Next, notice that

$\displaystyle\begin{split}\Delta_{g}u_{\epsilon}+\frac{C}{r^{2}}|\nabla
u_{\epsilon}|_{g}&=\Delta_{g}(\epsilon\psi)+\frac{C}{r^{2}}|\nabla(\epsilon\psi)|_{g}=\epsilon\Delta_{g}\psi+\epsilon\frac{C}{r^{2}}|\nabla\psi|_{g}\\\
&=\epsilon\left(\Delta_{g}\psi+\frac{C}{r^{2}}|\nabla\psi|_{g}\right)=-\frac{\epsilon}{r^{4}}\end{split}$
(3.2)

since $\psi$ satisfies (1.2) and we were able to move $\epsilon$ through the
norm since $\epsilon>0$. Then

$\displaystyle\begin{split}\mu_{\epsilon}&=u_{\epsilon}^{-4}\mu-4u_{\epsilon}^{-5}\Delta_{g}u_{\epsilon}=u_{\epsilon}^{-4}\mu+4u_{\epsilon}^{-5}\left(\frac{C}{r^{2}}|\nabla
u_{\epsilon}|_{g}+\frac{\epsilon}{r^{4}}\right)\\\
&>u_{\epsilon}^{-4}\mu_{\epsilon}+4u_{\epsilon}^{-5}|k_{ij}\nabla^{j}u_{\epsilon}|_{g}\geq
u_{\epsilon}^{-4}|J_{i}|_{{g}}+4u_{\epsilon}^{-5}|k_{ij}\nabla^{j}u_{\epsilon}|_{{g}}\\\
&\geq|(J_{\epsilon})_{i}|_{g_{\epsilon}}\end{split}$ (3.3)

where we used (2.6),(2.11), and (2.13). So indeed, the conformal initial data
set satisfies the strict dominant energy condition. Moreover, we can make the
estimate

$\displaystyle\mu_{\epsilon}-|J_{\epsilon}|_{g_{\epsilon}}\geq\frac{4\epsilon
u_{\epsilon}^{-5}}{r^{4}}>0$ (3.4)

and since for sufficiently small $\epsilon$ we have that $u_{\epsilon}$ is
close to $1$ we have

$\displaystyle\mu_{\epsilon}-|J_{\epsilon}|_{g_{\epsilon}}\geq\frac{\epsilon}{r^{4}}>0$
(3.5)

for all sufficiently small $\epsilon>0$.

The estimates (1.3) follow easily from the fact that $\psi\in
C^{2,\alpha}_{-1}({M})$ and the formula (2.2) for the ADM energy. Finally,
since $\psi$ satisfies the Neumann condition at the boundary, (2.20) shows
that

$\displaystyle(\theta_{\epsilon})_{\pm}=u_{\epsilon}^{-2}\theta_{\pm}$ (3.6)

on $\partial M$ so the sign of the null expansions is the same for any
$p\in\partial M$.

###  3.1. Solution of (1.2) on compact subdomains

The existence of a solution to (1.2) is established using similar methods to
Section 4.1 in [HirschKazarasKhuri] for the case of spacetime harmonic
functions. Let $S_{\mathcal{R}}$ denote a coordinate sphere of radius
$\mathcal{R}$ in $M$. Let $M_{\mathcal{R}}$ be the compact set defined by

$\displaystyle M_{\mathcal{R}}=\\{x\in M:r(x)\leq\mathcal{R}\\}$ (3.7)

which has boundary $\partial M_{\mathcal{R}}=\partial M\cup S_{\mathcal{R}}$.
Consider the mixed boundary value problem

$\displaystyle\begin{split}\Delta_{g}\psi_{\mathcal{R}}+\frac{C}{r^{2}}|\nabla\psi_{\mathcal{R}}|_{g}&=-\frac{1}{r^{4}}\quad\text{on}\quad
M_{\mathcal{R}}\\\
\left.\frac{\partial\psi_{\mathcal{R}}}{\partial\nu}\right|_{\partial
M}=0,\quad&\psi_{\mathcal{R}}|_{S_{\mathcal{R}}}=0\end{split}$ (3.8)

Recall we chose the coordinates (1.1) ensuring $r(x)\geq c_{0}>0$ for some
constant $c_{0}$ so that all the coefficients are smooth. We will solve this
subproblem using the Leray-Schauder fixed point theorem, see Theorem 11.3 in
[GilbargTrudinger]:

######  Theorem 3.1 (Leray-Schauder Fixed Point Theorem).

Let $\mathcal{B}$ be a Banach space and
$\mathcal{F}:\mathcal{B}\times[0,1]\rightarrow\mathcal{B}$ a compact mapping
with $\mathcal{F}(b,0)=0$ for all $b\in\mathcal{B}$. If there is a constant
$c$, such that any solution $(b,\sigma)\in\mathcal{B}\times[0,1]$ of
$b=\mathcal{F}(b,\sigma)$ satisfies the a priori inequality
$\|b\|_{\mathcal{B}}\leq c$, then there is a fixed point at $\sigma=1$. That
is, there exists $b_{1}\in\mathcal{B}$ with $b_{1}=\mathcal{F}(b_{1},1)$.

Define the space

$\displaystyle\widetilde{C}^{k,\alpha}({M}_{\mathcal{R}})=\\{f\in{C}^{k,\alpha}({M}_{\mathcal{R}}):\partial
f/\partial\nu|_{\partial M}=0,\quad f|_{S_{\mathcal{R}}}=0\\}$ (3.9)

consisting of functions satisfying the Neumann boundary condition on $\partial
M$ and the zero Dirichlet boundary condition on $S_{\mathcal{R}}$. Notice that
$\Delta^{-1}_{g}:C^{0,\alpha}\rightarrow\widetilde{C}^{2,\alpha}$ is an
isomorphism by applying the Fredholm alternative in the usual way. Next,
define the operator

$\displaystyle F(w)=-\frac{C}{r^{2}}|\nabla w|_{g}-\frac{1}{r^{4}}$ (3.10)

and notice if $w\in\widetilde{C}^{1,\alpha}(M_{\mathcal{R}})$ then $F(w)\in
C^{0,\alpha}(M_{\mathcal{R}})$. Letting
$\mathcal{B}=\widetilde{C}^{1,\alpha}(M_{\mathcal{R}})$ we define
$\mathcal{F}(w,\sigma):\mathcal{B}\times[0,1]\rightarrow\mathcal{B}$ to be the
composition of maps

$\displaystyle\widetilde{C}^{1,\alpha}(M_{\mathcal{R}})\xrightarrow{F}C^{0,\alpha}(M_{\mathcal{R}})\xrightarrow{\sigma\Delta_{g}^{-1}}\widetilde{C}^{2,\alpha}(M_{\mathcal{R}})\xrightarrow{i}\widetilde{C}^{1,\alpha}(M_{\mathcal{R}})$
(3.11)

which is a compact map since the first two maps are bounded and the last
inclusion is compact. Notice, $v=\mathcal{F}(w,\sigma)$ if $v$ solves the
problem

$\displaystyle\begin{split}&\Delta_{g}v=\sigma\left(-\frac{C}{r^{2}}|\nabla
w|_{g}-\frac{1}{r^{4}}\right)\\\ &\left.\frac{\partial
v}{\partial\nu}\right|_{\partial M}=0,\quad
v|_{S_{\mathcal{R}}}=0.\end{split}$ (3.12)

We denote the fixed points of the mapping by
$w_{\sigma}=\mathcal{F}(w_{\sigma},\sigma)$. Notice, $\mathcal{F}(w,0)=0$.
Next, we want to obtain an estimate of the form

$\displaystyle\|w_{\sigma}\|_{\widetilde{C}^{1,\alpha}(M_{\mathcal{R}})}\leq
c$ (3.13)

for all $0\leq\sigma\leq 1$. The key is to notice that we can always rewrite

$\displaystyle|\nabla w|_{g}=\frac{g_{ij}\nabla^{j}w}{|\nabla
w|_{g}}\nabla^{i}w$ (3.14)

where we can think of ${g_{ij}\nabla^{j}w}/{|\nabla w|_{g}}$ as a coefficient
in $L^{\infty}$. To be slightly more precise, this is undefined when $|\nabla
w|_{g}=0$ so we should define a coefficient depending on $\nabla w$ by

$\displaystyle b_{i}(\nabla
w,x)\coloneqq\begin{dcases}\frac{g_{ij}\nabla^{j}w}{|\nabla
w|_{g}}\quad&:\quad|\nabla w|_{g}\neq 0\\\ 0\quad&:\quad|\nabla
w|_{g}=0\end{dcases}.$ (3.15)

In this case indeed

$\displaystyle b_{i}(\nabla w,x)\nabla^{i}w=\begin{dcases}{|\nabla
w|_{g}}\quad&:\quad|\nabla w|_{g}\neq 0\\\ 0\quad&:\quad|\nabla
w|_{g}=0\end{dcases}\quad=|\nabla w|_{g}$ (3.16)

and $|b_{i}(\nabla w,x)|\leq c_{1}$ are uniformly bounded independent of
$\nabla w$. To improve readability, we write $b_{i}(\nabla w,x)=b_{i}(\nabla
w)$

Therefore, the fixed points of $\mathcal{F}$ satisfy the problem

$\displaystyle\begin{split}&\Delta_{g}w_{\sigma}=\sigma\left(-\frac{C}{r^{2}}b_{i}(\nabla
w_{\sigma})\nabla^{i}w_{\sigma}-\frac{1}{r^{4}}\right)\\\
&\left.\frac{\partial w_{\sigma}}{\partial\nu}\right|_{\partial M}=0,\quad
w_{\sigma}|_{S_{\mathcal{R}}}=0.\end{split}$ (3.17)

Therefore, if a solution exists, it can be viewed as a solution to the linear
problem

$\displaystyle\begin{split}&\Delta_{g}v=\sigma\left(-\frac{C}{r^{2}}b_{i}(\nabla
w_{\sigma})\nabla^{i}v-\frac{1}{r^{4}}\right)\\\ &\left.\frac{\partial
v}{\partial\nu}\right|_{\partial M}=0,\quad
v|_{S_{\mathcal{R}}}=0.\end{split}$ (3.18)

and this linear problem has uniformly bounded coefficients for all
$0\leq\sigma\leq 1$, independent of $w_{\sigma}$. Therefore, we can apply the
usual $L^{p}$ estimates for elliptic equations which give

$\displaystyle\|w_{\sigma}\|_{W^{2,p}(M_{\mathcal{R}})}\leq
C_{1}\left(\|r^{-4}\|_{{L^{p}(M_{\mathcal{R}})}}+\|w_{\sigma}\|_{L^{p}(M_{\mathcal{R}})}\right)$
(3.19)

for some elliptic constant $C_{1}$. Since the coefficient of the zero order
term in the linear PDE is zero, the maximum principle can be applied which
gives a $C^{0}$ and hence a $L^{p}$ bound for $w_{\sigma}$. Therefore we
obtain an estimate

$\displaystyle\|w_{\sigma}\|_{W^{2,p}(M_{\mathcal{R}})}\leq C_{2}$ (3.20)

for some constant $C_{2}$ independent of $\sigma$. For sufficiently large $p$,
the Sobolev embedding theorem gives ${W^{2,p}(M_{\mathcal{R}})}\hookrightarrow
C^{1,\alpha}(M_{\mathcal{R}})$ so we have the estimate

$\displaystyle\|w_{\sigma}\|_{C^{1,\alpha}(M_{\mathcal{R}})}\leq C_{3}$ (3.21)

independent of $\sigma$. We can now apply the Leray-Schauder fixed point
theorem for $\sigma=1$ to obtain a solution to (3.8) with at least
$\widetilde{C}^{1,\alpha}(M_{\mathcal{R}})$ regularity which we denote by
$w_{1}=\psi_{\mathcal{R}}$. But now we can look at it as as solution of

$\displaystyle\begin{split}\Delta_{g}\psi_{\mathcal{R}}&=f\quad\text{on}\quad
M_{\mathcal{R}}\\\
\left.\frac{\partial\psi_{\mathcal{R}}}{\partial\nu}\right|_{\partial
M}&=0,\quad\psi_{\mathcal{R}}|_{S_{\mathcal{R}}}=0\end{split}$ (3.22)

where

$\displaystyle
f=-\frac{C}{r^{2}}|\nabla\psi_{\mathcal{R}}|_{g}-\frac{1}{r^{4}}\in
C^{0,\alpha}(M_{\mathcal{R}})$ (3.23)

so therefore $\psi_{\mathcal{R}}\in C^{2,\alpha}(M_{\mathcal{R}})$.

Next, we need to show that we can take the limit
$\mathcal{R}\rightarrow\infty$. This can be done by constructing a pair of sub
and super solutions to (3.8) which are independent of $\mathcal{R}$. This is
done in the same way as in [HirschKazarasKhuri], except our case is even
simpler due to the form of our equation.

###  3.2. A uniform $C^{0}$ bound for $\psi_{\mathcal{R}}$

Let us remind ourselves what the definitions of sub and supersolutions are.
Given the equation

$\displaystyle Lu$ $\displaystyle=f(x)\quad\text{on}\quad\Omega$ (3.24)
$\displaystyle Bu$ $\displaystyle=h(x)\quad\text{on}\quad\partial\Omega$
(3.25)

where $L$ is a linear uniformly elliptic differential operator and
$Bu=A(x)\partial u/\partial\nu+B(x)u$ with either $A=0,B>0$ or $A>0,B\geq 0$
is a linear boundary operator, we say $\overline{u}$ is a supersolution if it
satisfies

$\displaystyle L\overline{u}$ $\displaystyle\leq
f(x)\quad\text{on}\quad\Omega$ (3.26) $\displaystyle B\overline{u}$
$\displaystyle\geq h(x)\quad\text{on}\quad\partial\Omega.$ (3.27)

Similarly we say $\underline{u}$ is a subsolution if it satisfies

$\displaystyle L\underline{u}$ $\displaystyle\geq
f(x)\quad\text{on}\quad\Omega$ (3.28) $\displaystyle B\underline{u}$
$\displaystyle\leq h(x)\quad\text{on}\quad\partial\Omega.$ (3.29)

If $\overline{u},\underline{u}\in H^{1}(\Omega)$ then we mean these
inequalities to hold in the weak sense.

Since we have shown $\psi_{\mathcal{R}}$ exists for each $\mathcal{R}$ we can
think of the family of linear problems with uniformly bounded coefficients

$\displaystyle\begin{split}L_{\mathcal{R}}(v)\coloneqq&\Delta_{g}v+\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}v=-\frac{1}{r^{4}}\quad\text{on}\quad
M_{\mathcal{R}}\\\ &\left.\frac{\partial v}{\partial\nu}\right|_{\partial
M}=0,\quad v|_{S_{\mathcal{R}}}=0.\end{split}$ (3.30)

and we would like to show there are sub and supersolutions to this problem
which are independent of $\mathcal{R}$. In fact, by the definition of a
subsolution we have $\underline{\psi}=0$ is a subsolution to (3.30) for any
$\mathcal{R}$. It remains to construct a supersolution $\overline{\psi}$.

Consider the function $\overline{w}(r)$

$\displaystyle\overline{w}(r)=\lambda
r^{-\beta},\quad\overline{w}^{\prime}(r)=-\lambda\beta
r^{-1-\beta},\quad\overline{w}^{\prime\prime}(r)=\lambda\beta(1+\beta)r^{-2-\beta}$
(3.31)

for $\lambda>0$ and $\beta\in(0,1)$. Using the fall-off conditions (2.1), it
can be shown that

$\displaystyle\Delta_{g}\overline{w}=-\lambda\beta(1-\beta)r^{-2-\beta}\left(1+O(r^{-1})\right),$
(3.32)

see equation (4.20) in [HirschKazarasKhuri]. It can also be verified that this
is correct up to the leading order term by calculating the Laplacian of
$\overline{w}$ with respect to the Euclidean metric.

Fix some $\beta\in(0,1)$. Then for any $\lambda>1$ there exists some $r_{0}$
such that

$\displaystyle\begin{split}&\Delta_{g}\overline{w}+\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}\overline{w}\leq-\frac{1}{r^{4}}\quad\text{for}\quad
r\geq r_{0}\end{split}$ (3.33)

independent of $\mathcal{R}$ (because for $\beta\in(0,1)$ the coefficient in
(3.32) is negative and the leading term decays as $r^{-2-\beta}$). Similar to
[HirschKazarasKhuri], we define the function

$\displaystyle\overline{\psi}=\begin{dcases}\widetilde{\psi}\coloneqq\widehat{\psi}_{r_{0}}+\lambda
r_{0}^{-\beta}\quad&:\quad r\leq r_{0}\\\ \overline{w}=\lambda
r^{-\beta}\quad&:\quad r>r_{0}\end{dcases}$ (3.34)

where $\widehat{\psi}_{r_{0}}$ is the solution of (3.8) on $M_{r_{0}}$ but
with $-2r^{-4}$ instead of $-r^{-4}$ on the right hand side of the equation.
Notice, $\lambda r_{0}^{-\beta}$ is just a constant. We see that
$\overline{\psi}$ is $C^{2,\alpha}$ on $M_{r_{0}}$, smooth on $M\setminus
M_{r_{0}}$ and by construction continuous on $S_{r_{0}}$ for any $\lambda>1$.
It is therefore easy to show $\overline{\psi}\in H^{1}(M_{\mathcal{R}})$ for
any $\mathcal{R}$.

Let $\partial_{r}f$ denote the outward normal derivative of $f$ on
$S_{r_{0}}$. Due to (3.31), for a sufficiently large $r_{0}$ we can take a
sufficiently large $\lambda>1$ such that

$\displaystyle\partial_{r}\widetilde{\psi}>\partial_{r}\overline{w}\quad\text{on}\quad
S_{r_{0}}.$ (3.35)

We can now show that $\overline{\psi}$ is a supersolution. Notice,

$\displaystyle b_{i}(\nabla w)\nabla^{i}v\leq|\nabla v|_{g}$ (3.36)

for any $w,v$. Therefore

$\displaystyle\begin{split}L_{\mathcal{R}}(\widetilde{\psi})&=\Delta_{g}\widetilde{\psi}+\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}\widetilde{\psi}\leq\Delta_{g}\widetilde{\psi}+\frac{C}{r^{2}}|\nabla\widetilde{\psi}|_{g}=\Delta_{g}\widehat{\psi}_{r_{0}}+\frac{C}{r^{2}}|\nabla{\widehat{\psi}_{r_{0}}}|_{g}=-\frac{2}{r^{4}}\\\
&<-\frac{1}{r^{4}}\end{split}$ (3.37)

on $M_{r_{0}}\cap M_{\mathcal{R}}$. Therfore $\overline{\psi}$ is a
supersolution on $M_{r_{0}}\cap M_{\mathcal{R}}$ and $M_{\mathcal{R}}\setminus
M_{r_{0}}$ separately. The problem is that the function is not differentiable
across $S_{r_{0}}$. However, we can show that it is a weak supersolution which
is enough.

For our problem, $\overline{\psi}$ will be a supersolution if for a non-
negative test function $\phi\in C_{c}^{\infty}(M_{\mathcal{R}})$ we have

$\displaystyle\int_{M_{\mathcal{R}}}\left(-\langle\nabla\phi,\nabla\overline{\psi}\rangle_{g}+\phi\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}\overline{\psi}+\phi\frac{1}{r^{4}}\right)dV\leq
0.$ (3.38)

So let $\phi\in C_{c}^{\infty}(M_{\mathcal{R}})$ be a non-negative test
function. By (3.37) since
$L_{\mathcal{R}}(\widetilde{\psi})\leq-r^{-4}=L_{\mathcal{R}}(\psi_{\mathcal{R}})$
we have

$\displaystyle\begin{split}0&\geq\int_{M_{\mathcal{R}}\cap
M_{r_{0}}}\phi(L_{\mathcal{R}}(\widetilde{\psi})+r^{-4})dV=\int_{M_{\mathcal{R}}\cap
M_{r_{0}}}\phi(L_{\mathcal{R}}(\overline{\psi})+r^{-4})dV\\\
&=\int_{M_{\mathcal{R}}\cap
M_{r_{0}}}\phi\left(\Delta_{g}\overline{\psi}+\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}\overline{\psi}+\frac{1}{r^{4}}\right)dV\\\
&=\int_{M_{\mathcal{R}}\cap
M_{r_{0}}}\left(-\langle\nabla\phi,\overline{\psi}\rangle_{g}+\phi\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}\overline{\psi}+\phi\frac{1}{r^{4}}\right)dV+\int_{S_{r_{0}}}\phi\partial_{r}\overline{\psi}dA\\\
&=\int_{M_{\mathcal{R}}\cap
M_{r_{0}}}\left(-\langle\nabla\phi,\overline{\psi}\rangle_{g}+\phi\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}\overline{\psi}+\phi\frac{1}{r^{4}}\right)dV+\int_{S_{r_{0}}}\phi\partial_{r}\widetilde{\psi}dA\end{split}$
(3.39)

where we integrated by parts. Notice there is no boundary integral term on
$\partial M$ since $\phi\in C_{c}^{\infty}(M_{\mathcal{R}})$ and the last term
has a positive sign because the outward unit normal to $M_{\mathcal{R}}\cap
M_{r_{0}}$ at $S_{r_{0}}$ points to infinity.

Similarly using (3.33) we have

$\displaystyle\begin{split}0&\geq\int_{M_{\mathcal{R}}\setminus
M_{r_{0}}}\phi(L_{\mathcal{R}}(\overline{w})+r^{-4})dV=\int_{M_{\mathcal{R}}\setminus
M_{r_{0}}}\phi(L_{\mathcal{R}}(\overline{\psi})+r^{-4})dV\\\
&=\int_{M_{\mathcal{R}}\setminus
M_{r_{0}}}\phi\left(\Delta_{g}\overline{\psi}+\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}\overline{\psi}+\frac{1}{r^{4}}\right)dV\\\
&=\int_{M_{\mathcal{R}}\setminus
M_{r_{0}}}\left(-\langle\nabla\phi,\overline{\psi}\rangle_{g}+\phi\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}\overline{\psi}+\phi\frac{1}{r^{4}}\right)dV-\int_{S_{r_{0}}}\phi\partial_{r}\overline{\psi}dA\\\
&=\int_{M_{\mathcal{R}}\setminus
M_{r_{0}}}\left(-\langle\nabla\phi,\overline{\psi}\rangle_{g}+\phi\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}\overline{\psi}+\phi\frac{1}{r^{4}}\right)dV-\int_{S_{r_{0}}}\phi\partial_{r}\overline{w}dA\end{split}$
(3.40)

where again there is no surface integral on $S_{M_{\mathcal{R}}}$ because
$\phi\in C_{c}^{\infty}(\mathcal{R})$ and the last term has a minus sign
because the outward unit normal to $M_{\mathcal{R}}\setminus M_{r_{0}}$ on
$S_{r_{0}}$ points away from infinity.

Adding these two inequalities together then gives

$\displaystyle
0\geq\int_{M_{\mathcal{R}}}\left(-\langle\nabla\phi,\nabla\left(\overline{\psi}\right)\rangle_{g}+\phi\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}(\overline{\psi})+\phi\frac{1}{r^{4}}\right)dV+\int_{S_{r_{0}}}\phi\partial_{r}(\widetilde{\psi}-\overline{w})dA$
(3.41)

or

$\displaystyle\int_{S_{r_{0}}}\phi\partial_{r}(\overline{w}-\widetilde{\psi})dA\geq\int_{M_{\mathcal{R}}}\left(-\langle\nabla\phi,\nabla\left(\overline{\psi}\right)\rangle_{g}+\phi\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}(\overline{\psi})+\phi\frac{1}{r^{4}}\right)dV.$
(3.42)

But, by (3.35) the left hand side is non-positive so therefore

$\displaystyle
0\geq\int_{M_{\mathcal{R}}}\left(-\langle\nabla\phi,\nabla\left(\overline{\psi}\right)\rangle_{g}+\phi\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}(\overline{\psi})+\phi\frac{1}{r^{4}}\right)dV$
(3.43)

and so indeed $\overline{\psi}$ is a weak supersolution to (3.30) for any
$\mathcal{R}$.

Since $\widehat{\psi}_{r_{0}}$ in the definition of $\overline{\psi}$ is a
fixed function, we can take a sufficiently large $\lambda$ such that
$\overline{\psi}>0$ everywhere. We’d like to show
$\overline{\psi}\geq\psi_{\mathcal{R}}$ for all $\mathcal{R}$ to obtain a
uniform $C^{0}$ bound for all the solutions.

Let $V_{\mathcal{R}}=\overline{\psi}-\psi_{R}$ on $M_{\mathcal{R}}$ so that
$L(V_{\mathcal{R}})\leq 0$ on $M_{\mathcal{R}}$. In fact we have

$\displaystyle
L_{\mathcal{R}}(V_{\mathcal{R}})=L_{\mathcal{R}}(\overline{\psi})-L_{\mathcal{R}}(\psi_{\mathcal{R}})\leq-\frac{2}{r^{4}}+\frac{1}{r^{4}}=-\frac{1}{r^{4}}<0\quad\text{on}\quad
M_{\mathcal{R}}\cap M_{r_{0}}.$ (3.44)

Since $V_{\mathcal{R}}\in H^{1}(M_{\mathcal{R}})\cap C^{0}(M_{\mathcal{R}})$
the weak maximum principle (see Theorem 8.1 in [GilbargTrudinger]) gives

$\displaystyle\inf_{M_{\mathcal{R}}}V_{\mathcal{R}}=\inf_{\partial
M_{\mathcal{R}}}V_{\mathcal{R}}.$ (3.45)

If $\inf_{\partial M_{\mathcal{R}}}V_{\mathcal{R}}\geq 0$ we’re done. So
suppose $V_{\mathcal{R}}<0$ somewhere, so that the infimum occurs necessarily
somewhere on $\partial M_{\mathcal{R}}$. Now, $\partial
M_{\mathcal{R}}=\partial M\cup S_{\mathcal{R}}$. On $S_{\mathcal{R}}$ we have
$V_{\mathcal{R}}>0$ because $\overline{\psi}>0$ and $\psi_{\mathcal{R}}=0$ on
$S_{\mathcal{R}}$. Therefore, let us suppose that the infimum is reached at
some point $x_{0}$ on $\partial M$.

Now take $\overline{r}=\min\\{r_{0},\mathcal{R}\\}$. Consider the set
$\Lambda=M_{\overline{r}}$ and its interior $\Lambda^{o}$. We have
$V_{\mathcal{R}}\in C^{2,\alpha}(\Lambda)$ and the infimum on this set is
still reached at the same $x_{0}\in\partial M$. We can now apply Lemma 3.4
from [GilbargTrudinger] which gives that if
$V_{\mathcal{R}}(x_{0})<V_{\mathcal{R}}(x)$ for all $x\in\Lambda^{o}$ then
necessarily

$\displaystyle\frac{\partial V_{\mathcal{R}}}{\partial\nu}(x_{0})<0.$ (3.46)

But this is a contradiction, since $x_{0}\in\partial M$ and our functions
satisfy the zero Neumann condition there. (We remark that this boundary point
lemma is proven for the Euclidean normal derivative in the given choice of
coordinates, but it can be easily generalized to apply to the normal
derivative with respect to the metric. Another way to see this is to take
geodesic normal coordinates at $x_{0}$ so that the two normal derivatives
match at that point in these coordinates.)

Therefore, there must be some point $x\in\Lambda^{o}$ such that
$V_{\mathcal{R}}(x_{0})=V_{\mathcal{R}}(x)$. But then we can apply the strong
minimum principle to conclude $V_{\mathcal{R}}$ is constant on $\Lambda\subset
M_{\mathcal{R}}\cap M_{r_{0}}$. But this then contradicts (3.44).

Thus, we concldue that the infimum of $V_{\mathcal{R}}$ is not reached on
$\partial M$. Therefore, it is reached on $S_{\mathcal{R}}$ where
$V_{\mathcal{R}}$ is positive. Therefore

$\displaystyle\inf_{M_{\mathcal{R}}}(\overline{\psi}-\psi_{\mathcal{R}})=\inf_{M_{\mathcal{R}}}V_{\mathcal{R}}=\inf_{\partial
M_{\mathcal{R}}}V_{\mathcal{R}}\geq 0$ (3.47)

so that $\psi_{\mathcal{R}}\leq\overline{\psi}$. It is also elementary to show
$0\leq\psi_{\mathcal{R}}$ giving us the uniform $C^{0}$ estimate

$\displaystyle 0\leq\psi_{\mathcal{R}}\leq\overline{\psi}.$ (3.48)

###  3.3. Global existence

We can now use the uniform estimate (3.48) to construct a convergent
subsequence to obtain a global solution of (1.2). Take a countable cover of
$M$ consisting of smooth precompact open sets $\\{\Omega_{n}\\}$ (for example
a countable collection of balls of Euclidean radius $1$ in our asymptotically
flat coordinates intersected with $M$ will do).

We can think of the functions $\psi_{\mathcal{R}}$ as solutions of the family
of elliptic operators with uniformly bounded coefficients as in
$\eqref{MixedProblemSubSuper}$ but now on each $\Omega_{n}$. For each
$\Omega_{n}$ we have $\Omega_{n}\subset M_{\mathcal{R}}$ for all sufficiently
large $\mathcal{R}$ due to the precompactness of the elements of the cover.

Therefore we can apply the usual $L^{p}$ estimates once again to obtain

$\displaystyle\|\psi_{\mathcal{R}}\|_{W^{2,p}(\Omega_{n})}\leq
C_{n}\left(\|r^{-4}\|_{{L^{p}(\Omega_{n})}}+\|\psi_{\mathcal{R}}\|_{L^{p}(\Omega_{n})}\right)$
(3.49)

where the elliptic constant $C_{n}$ depends on $\Omega_{n}$ but not on
$\mathcal{R}$ due to the uniform boundedness of the coefficients in (3.30).
Combined with (3.48) this then gives

$\displaystyle\|\psi_{\mathcal{R}}\|_{W^{2,p}(\Omega_{n})}\leq
C_{n}\left(\|r^{-4}\|_{{L^{p}(\Omega_{n})}}+\|\overline{\psi}\|_{L^{p}(\Omega_{n})}\right)\leq
C^{\prime}_{n}$ (3.50)

for some constant $C^{\prime}_{n}$ independent of $\mathcal{R}$.

Choosing suffieiently large $p$ and applying the Sobolev embedding theorem
then gives a uniform estimate

$\displaystyle\|\psi_{\mathcal{R}}\|_{C^{1,\gamma}(\overline{\Omega}_{n})}\leq
C_{n}^{\prime\prime}$ (3.51)

for any $0<\gamma<1$ for some constant $C_{n}^{\prime\prime}$ independent of
$\mathcal{R}$. Using this estimate and looking at the $\psi_{\mathcal{R}}$ as
solutions of

$\displaystyle\begin{split}&\Delta_{g}\psi_{\mathcal{R}}=f_{\mathcal{R}}\quad\text{on}\quad\Omega_{n}\\\
f_{\mathcal{R}}&=-\frac{C}{r^{2}}|\nabla\psi_{\mathcal{R}}|_{g}-\frac{1}{r^{4}}\in
C^{0,\gamma}(\overline{\Omega}_{n})\end{split}$ (3.52)

we have that the $f_{\mathcal{R}}$ are uniformly bounded in
$C^{0,\gamma}(\overline{\Omega}_{n})$. Combining this with the uniform
estimate $\eqref{C0Estimate}$ and the usual Schauder estimates gives the
uniform estimate

$\displaystyle\|\psi_{\mathcal{R}}\|_{C^{2,\gamma}(\overline{\Omega}_{n})}\leq
C_{n}^{\prime\prime\prime}$ (3.53)

Using the fact that the embedding $C^{2,\gamma}\hookrightarrow C^{2,\alpha}$
is compact for $0<\alpha<\gamma<1$, we can take a sequence
$\psi_{\mathcal{R}_{i}}$ with $\mathcal{R}_{i}\rightarrow\infty$ which on any
$\Omega_{n}$ has a convergent subsequence converging in
$C^{2,\alpha}(\overline{\Omega}_{n})$. Using the usual diagonal argument, we
can then extract a subsequence converging on all $\Omega_{n}$ and thus obtain
a subsequence converging on compact subsets of $M$ to a $C^{2,\alpha}(M)$
solution $\psi$ of (1.2) satisfying

$\displaystyle 0\leq\psi\leq\lambda r^{-\beta}\quad\text{for}\quad r\geq
r_{0}.$ (3.54)

###  3.4. Estimates at infinity

Next, we want to show that the resulting solution $\psi$ is in fact in the
weighted Hölder space $C^{2,\alpha}_{-1}(M)$. In [HirschKazarasKhuri], the
authors refer to the methods of Proposition 3 in [SchoenYau2] in order to
obtain their desired fall-off conditions, but a careful look at that
proposition only shows the solution to be in $C^{2,\alpha}_{-\beta}$ for
$\beta\in(0,1)$ which is not quite good enough for us. Thus, we make use of
weighted elliptic estimates from [Chaljub1979] and [Weinstein_2005].

Consider the region

$\displaystyle\Omega=\\{x\in M:r(x)\geq\overline{r}_{0}>r_{0}\\}.$ (3.55)

Take some point $x\in\Omega$ with $r(x)=r_{x}$ such that the ball of Euclidean
radius $1$ centered at $x$ is contained in $\Omega$ so
$B=B_{1}(x)\subset\Omega$. Notice, if we take $\overline{r}_{0}$ sufficiently
large, then we have

$\displaystyle\frac{1}{r_{x}-1}\leq\frac{2}{r_{x}}$ (3.56)

so that on $B$ we have

$\displaystyle\frac{1}{r(y)}\leq\frac{2}{r_{x}}$ (3.57)

for all $y\in B$.

Now on each such ball, we can look at $\psi$ as a solution to the problem

$\displaystyle\begin{split}L(v)\coloneqq&\Delta_{g}v+\frac{C}{r^{2}}b_{i}(\nabla\psi_{\mathcal{R}})\nabla^{i}v=-\frac{1}{r^{4}}\quad\text{on}\quad
B\end{split}$ (3.58)

which as before allows us to use the usual $L^{p}$ estimates giving

$\displaystyle\|\psi\|_{W^{2,p}(B)}\leq
C_{e}\left(\|r^{-4}\|_{{L^{p}(B)}}+\|\psi\|_{L^{p}(B)}\right)$ (3.59)

where the elliptic constant $C_{e}$ can be chosen so that it does not depend
on the center of the ball due to the uniform boundedness of the coefficients
(which is also a consequence of the asymptotic flatness). Since the ball has
radius $1$ and is in $\Omega$ then we can estimate

$\displaystyle\left(\int_{B}|r^{-4}|^{p}dV\right)^{1/p}\leq\left(\int_{B}|(2/r_{x})^{4}|^{p}dV\right)^{1/p}=16\left(\frac{4}{3}\pi\right)^{1/p}r_{x}^{-4}$
$\displaystyle\left(\int_{B}|\psi|^{p}dV\right)^{1/p}\leq\left(\int_{B}|\lambda
r^{-\beta}|^{p}dV\right)^{1/p}\leq\left(\int_{B}|\lambda
2^{\beta}r_{x}^{-\beta}|^{p}dV\right)^{1/p}=\lambda
2^{\beta}\left(\frac{4}{3}\pi\right)^{1/p}r_{x}^{-\beta}.$

Putting these estimates together we then have

$\displaystyle\|\psi\|_{W^{2,p}(B)}\leq C_{e}^{\prime}r_{x}^{-\beta}$ (3.60)

for some constant $C_{e}^{\prime}$ depending on $\lambda,\beta,p$ but not the
center of the ball $x$. For sufficiently large $p$ we can now apply the
Sobolev embedding theorem, where the Sobolev constant does not depend on $x$
since we are looking at Euclidean balls of radius $1$. Therefore we get

$\displaystyle\|\psi\|_{C^{1,\gamma}(\overline{B})}\leq
C_{e}^{\prime\prime}r_{x}^{-\beta}$ (3.61)

where $\gamma=1-3/p$ and $C_{e}^{\prime\prime}$ doesn’t depend on the center
of the ball. Notice, this means $\partial_{i}\psi=O(r^{-\beta})$.

Now, let us look at

$\displaystyle f=-\frac{C}{r^{2}}|\nabla\psi|_{g}-\frac{1}{r^{4}}\in
C^{0,\gamma}(\overline{B})$ (3.62)

which satisfies

$\displaystyle\|f\|_{C^{0,\gamma}(\overline{B})}\leq
C_{e}^{\prime\prime\prime}r_{x}^{-\beta-2}$ (3.63)

for some constant $C_{e}^{\prime\prime\prime}$ independent of $x\in\Omega$ by
(3.61).

Now, take some $0<\alpha^{\prime}<\gamma$ such that
$\alpha^{\prime}+\beta^{\prime}<\beta$. Now, take any $y\in B_{1/4}(x)$. Let
$\widetilde{r}(x,y)=\min\\{r(x),r(y)\\}<r_{x}$. We have

$\displaystyle\begin{split}\sup_{y\in
B_{1/4}(x)}(\widetilde{r}^{2+\beta^{\prime}+\alpha^{\prime}}(x,y))\frac{|f(x)-f(y)|}{|x-y|^{\alpha^{\prime}}}&\leq
r_{x}^{2+\beta^{\prime}+\alpha^{\prime}}\left(\sup_{y\in
B_{1/4}(x)}\frac{|f(x)-f(y)|}{|x-y|^{\alpha^{\prime}}}\right)\\\ &\leq
r_{x}^{2+\beta^{\prime}+\alpha^{\prime}}\left(\sup_{y\in
B_{1/4}(x)}\frac{|f(x)-f(y)|}{|x-y|^{\gamma}}\right)\\\ &\leq
r_{x}^{2+\beta^{\prime}+\alpha^{\prime}}\left(C_{e}^{\prime\prime\prime}r_{x}^{-\beta-2}\right)<\infty\end{split}$
(3.64)

Therefore by definition of weighted Hölder spaces (again, we refer the reader
to [EichmairHuangLeeSchoen2016, Chaljub1979, Weinstein_2005] for the
definitions) we have $f\in
C^{0,\alpha^{\prime}}_{-\beta^{\prime}-2}(\overline{\Omega})$ and so in fact

$\displaystyle f\in C^{0,\alpha^{\prime}}_{-\beta^{\prime}-2}(M).$ (3.65)

Looking at $\psi$ as a solution of

$\displaystyle\begin{split}&\Delta_{g}v=f\quad\text{on}\quad M\\\
&\left.\frac{\partial v}{\partial\nu}\right|_{\partial M}=0\end{split}$ (3.66)

we can now use the results of [Chaljub1979] to conclude that there exists a
unique solution to (3.66) in
$\widetilde{C}^{2,\alpha^{\prime}}_{-\beta^{\prime}}(M)$, which therefore is
the same $\psi$ we found earlier.

######  Remark 1.

Technically, the results of [Chaljub1979] are proven for the case of complete
asymptotically flat manifolds without boundary, but as mentioned in the
footnote on page 11 of that reference, the results can be generalized to the
case of manifolds with boundary. In particular, this is easy to see in our
case since we can apply the Fredholm alternative to (3.66).

Since we now have $\partial_{i}\psi=O(r^{-\beta^{\prime}-1})$ this improves
the asymptotics of $f$ giving

$\displaystyle f\in C^{0,\alpha^{\prime}}_{-2-\beta}$ (3.67)

for any $\beta\in(0,1)$, yielding that
$\psi\in\widetilde{C}^{2,\alpha^{\prime}}_{-\beta}(M)$. Repeating this
argument one final time again improves the asymptotics on $f$ and thus on
$\psi$ giving

$\displaystyle\psi\in\widetilde{C}^{2,\alpha}_{-\beta}(M)$ (3.68)

for any $\alpha,\beta\in(0,1)$.

Now, notice if we take $\beta=-2/3$ then

$\displaystyle f=-\frac{C}{r^{2}}|\nabla\psi|_{g}-\frac{1}{r^{4}}\in
C^{0,\alpha}_{-3}\cap L^{1}(M).$ (3.69)

We can now use the elliptic estimate from Lemma 2 of [Weinstein_2005] which
gives

$\displaystyle\|\psi\|_{C^{2,\alpha}_{-1}(M)}\leq
C\left(\|f\|_{C^{0,\alpha}_{-3}\cap
L^{1}(M)}+\|\psi\|_{C^{2,\alpha}_{-2/3}(M)}\right)$ (3.70)

for some constant $C$. Therefore indeed $\psi\in C^{2,\alpha}_{-1}(M)$ as
desired.

##  4\. Possible Further Applications

One possible application of the methods of this paper is to the deformation of
charged initial data sets. One can consider an initial data set $(M,g,k,E,B)$
where $E$ and $B$ are the electric and magnetic fields. Then the charged
dominant energy condition takes the form

$\displaystyle\mu_{EM}\geq|J_{EM}|_{g}$ (4.1)

where $\mu_{EM}$ and $J_{EM}$ denote the energy and momentum densities
respectively with the contributions from the electromagnetic fields subtracted
off. As far as we have seen, no general result about deforming such initial
data sets to the strict charged dominant energy condition have been published,
and the results of the current paper should readily be able to handle this
case.

Some other possible applications are to the Penrose inequality, which is one
of the most important open problem in mathematical general relativity
[Penrose], which relates the total mass of a spacetime to the area of the
black holes present in the spacetime. The mathematically precise formulation
of the conjecture for asymptotically flat initial data sets is believed to be

$\displaystyle E_{ADM}=\sqrt{\frac{A}{16\pi}}$ (4.2)

where $A$ is the total area of the outermost area enclosure of future/past
apparent horizons present in the initial data sets. The inequality has been
proven in the so-called Riemannian case where $k=0$ and the dominant energy
condition reduces to non-negative scalar curvature, and the outermost apparent
horizon reduces to an outermost minimal surface [Bray, HuiskenIlmanen].
Various generalizations to include charge and angular momentum have been
proposed and proven under certain conditions on $k$ [KhuriWeinsteinYamada,
KhuriWeinsteinYamada1, JaraczKhuri, JaraczCylindricalPenrose,
KhuriSokolowskyWeinstein], but even the simplest version (4.2) for general $k$
remains open. While some general approaches have been proposed [BrayKhuriPDE],
it has turned out some of these approaches do not work, as for example the
proposed Jang/zero-divergence system [Jaracz_2023_ZeroDivergence].

If the Penrose conjecture turns out to be true, then the dominant energy
condition is a key ingredient. In fact the heuristic physics argument which
leads to the Penrose inequality relies on the black hole area theorem which
holds under the weaker assumption of either the weak energy condition or the
null energy condition. The dominant energy condition implies the weak energy
condition which in turn implies the null energy condition, but not conversely.
However, there exist counterexamples to both the positive mass theorem and the
Penrose inequality under the assumption of the weak energy condition
[JaraczWEC]. It remains unclear what causes this breakdown between the
physical and mathematical arguments.

The most powerful approach to the Penrose inequality was the conformal
approach of Bray [Bray] where the flow was defined by

$\displaystyle\begin{split}u_{t}&=1+\int_{0}^{t}v_{s}ds\\\
\Delta_{g}v_{t}&=0\quad\text{on}\quad M_{t},\\\ v_{t}\equiv
0\quad\text{on}\quad M\setminus M_{t},\quad v_{t}|_{\partial
M_{t}}&=0,\quad\lim_{|x|\rightarrow\infty}v_{t}(x)=-e^{-t}.\end{split}$ (4.3)

Here, $M_{t}$ is the region outside the outermost minimal surface in
$(M,g_{t})$ where $g_{t}=u_{t}^{4}g$. This conformal flow preserved the non-
negativity of scalar curvature.

Similarly, looking at (2.12) we see that if we have a family of conformal
factors which satisfy

$\displaystyle\Delta_{g}u_{t}=-|k_{ij}\nabla^{j}u_{t}|_{g}$ (4.4)

then this family of metrics will preserve the dominant energy condition. In
fact, we could have studied this type of equation instead, but using the
estimate (2.13) proved to be more convenient in the present work. In addition,
if we assume

$\displaystyle u_{t}$ $\displaystyle=1+\int_{0}^{t}v_{s}ds$ (4.5)

so that $u_{t}^{\prime}=v_{t}$ we see that (4.4) is true for $t=0$. Then
letting $k_{ij}\nabla^{j}u_{t}=Y_{i}$ and differentiating (4.4) with respect
to $t$ gives

$\displaystyle\Delta_{g}v_{t}=-\frac{Y^{i}k_{ij}}{|Y_{i}|}\nabla^{j}v_{t}.$
(4.6)

In the special case where $k\equiv 0$ this reduces to (4.3). It is unclear if
such a flow could be used to tackle the general case of the Penrose
inequality, but it looks like a promising possibility.

## References

J.S. Jaracz, Department of Mathematics, Texas State University, San Marcos, TX
78666

E-mail address<EMAIL_ADDRESS>


chopped. Chopping helps to suppress the sky background, which gets
increasingly dominant in the $M$ band. Thus, in $M$ band it is mandatory to
use chopped data, while in $L$ band both chopped and non-chopped data may be
used. During chopping the frequent switching between the object and sky
positions might cause to lose track of the adaptive optics guiding or the
fringe tracking, which results in bad quality data. Thus, when both non-
chopped and chopped $L$ band data-sets were available, we carefully inspected
both products, and selected the one which we deemed more reliable (the one
with smaller error bars).

##  Appendix C Visibility calculation

The visibility function is the Fourier-transform of the object’s image on the
sky. IR interferometers sample the visibility function at distinct spatial
frequencies, set by the baseline length between the telescope pairs. The
visibility function is a complex quantity, and in IR interferometry it is
usually represented by its absolute value and phase. The absolute value of the
visibility can be expressed either as a unitless quantity, or in flux units.
In the latter case it is known as the correlated flux. It is the correlated
flux which is directly measured by the interferometer. In polar coordinates
the Fourier-transform becomes the Hankel-transform. For an object with
circular symmetry, the correlated flux can be expressed from the radial
surface brightness profile ($I_{\nu}\left(r\right)$) as follows:

$F_{\text{corr},\nu}\left(B_{\mathrm{p}}\right)=\int_{R_{\text{in}}}^{R_{\text{out}}}2\pi
rI_{\nu}\left(r\right)J_{0}\left(2\pi
rB_{\mathrm{p}}/\lambda\right)\mathrm{d}r,$ (8)

where $B_{\text{p}}$ is the projected baseline length, $J_{0}$ is the Bessel
function of the first kind of order zero, and $R_{\text{in}}$,
$R_{\text{out}}$ are the inner and outer radii of the object, respectively.
The $B_{\text{p}}/\lambda$ is the spatial frequency. The argument of $J_{0}$
is in radians. In our multi-zone disk model, $F_{\text{corr},\nu,i}$ is the
correlated flux of the i-th zone, with the integration limits being the zone
boundaries $R_{\text{in}}=R_{\text{i}}$ and $R_{\text{out}}=R_{\text{i+1}}$.
The formula above assumes a circularly symmetric object image. If the image is
elliptic, such as an inclined thin disk, the visibility function also depends
on the position angle of the baselines ($\theta_{B}$), not just on their
lengths. In that case we can transform the projected baseline lengths so that
we remove the effect of the inclination:

$B_{\text{eff}}=B_{p}\sqrt{\cos^{2}\left(\Delta\theta\right)+\cos^{2}\left(i\right)\sin^{2}\left(\Delta\theta\right)},$
(9)

where $\Delta\theta=\theta_{B}-\theta$ is the difference between the position
angle of the baseline and the position angle of the major axis of the object
image ($\theta$). Using $B_{\text{eff}}$ removes the $\theta_{B}$ dependence
in the visibility calculation, and we can apply Eq. 8, as if the object image
were circular.

In our model the central star is represented as a point source, thus its
correlated flux equals its (single-dish) flux. The correlated fluxes of the
disk radial zones ($F_{\text{corr},\nu,i}$), and that of the central star
($F_{\nu,\text{star}}$) are simply added together:

$F_{\text{corr},\nu}\left(B_{\mathrm{eff}}\right)=\sum_{i=1}^{M}F_{\text{corr},\nu,i}\left(B_{\mathrm{eff}}\right)+F_{\nu,\text{star}}.$
(10)

Finally, we calculate the unitless visibility ($V_{\nu}$), by normalizing the
correlated flux with the single-dish flux ($F_{\nu}$):

$V_{\nu}\left(B_{\mathrm{eff}}\right)=F_{\text{corr},\nu}\left(B_{\mathrm{eff}}\right)/F_{\nu}.$
(11)

##  Appendix D Opacity curves

Figure 15: $7-35\ \mu$m opacity curves of pyroxene, olivine, and various iron
oxides, compared with the normalized Spitzer spectrum of HD 144432 (brown
curve, AOR: 3587072). The Spitzer spectrum was normalized by dividing the data
with the source function of our best-fit model. Optical constants of iron
oxides were downloaded from the Jena Database of Optical Constants for Cosmic
Dust (https://www.astro.uni-jena.de/Laboratory/OCDB/mgfeoxides.html). The
original source for wüstite is Henning et al. (1995). Optical data for
hematite and magnetite come from Amaury H.M.J. Triaud (unpublished). The
opacity curves were generated with optool, using the distribution of hollow
spheres (DHS) scattering theory. The upper boundary of the hollow sphere
distribution ($f_{\mathrm{max}}$) is $0.7$. All curves are arbitrarily scaled
and shifted to aid visual comparison.  Figure 16: $N$ band opacity curves of
the silicate minerals used in this work. Red color indicates the subset used
in the final dust modeling.  Figure 17: Opacity curves in the full wavelength
range of the data. Shown only the subset used in the final dust modeling,
including iron and carbon as well.

##  Appendix E Two-step optimization procedure

Our approach for finding the best-fit model is detailed in the following.
First, we constrain the silicate composition and the main structural disk
parameters by fitting only the $N$ band MATISSE data (six correlated spectra
and the single-dish spectrum between $8$ and $13\ \mu$m). For this purpose we
use our own idl software, based on the `specfit` dust compositional fitting
tool by Juhász et al. (2009). In our code we implemented our multi-zone
temperature gradient model, and added the functionality to model
interferometric data, while we retained the original fitting procedure of
`specfit`. The optimization of the disk structural parameters is done by a
genetic optimization algorithm named `pikaia` (Charbonneau 1995). In each
optimization step of the genetic algorithm, a spectral decomposition is
employed by a bound constrained linear least-squares fitting algorithm (Lawson
& Hanson 1974), in order to evaluate the dust mass fractions. The latter step
requires a change in our model when calculating the surface brightness profile
in Eq. 4:

$I_{\nu,i}\left(r\right)=\tau_{\nu,i}\left(r\right)B_{\nu}\left(T\left(r\right)\right)=\sum_{j=1}^{N}{c_{i,j}\Sigma_{i}\left(r\right)\kappa_{\nu,j}}B_{\nu}\left(T\left(r\right)\right).$
(12)

Here we assume the disk to be optically thin ($\tau\ll 1$) at our observing
wavelengths. Thus, we can calculate the total emission by simple linear
addition of the contributions from the individual dust components. This is
needed by the spectral decomposition algorithm in the specfit code. The final
result of this modeling is the set of fit structural parameters and dust mass
fractions.

We expect that the N-band data alone cannot differentiate between the dust
species with featureless opacity curves (iron and carbon). This is because the
opacities of those species look very similar in that spectral region. Thus, we
do not include either carbon or iron in our `specfit` runs. Instead we include
a hidden dust component with a power law opacity curve to represent the dust
responsible for the featureless continuum. The exponent of the power law is a
fit parameter, while we fix the opacity191919This represents a typical order-
of-magnitude value for submicron-sized grains. at $8\ \mu$m at $500\
\text{cm}^{2}/\text{g}$.

In the second modeling step we aim to constrain the disk structure overall.
The fitting procedure employed in the first step is very capable at
constraining the dust composition, but it can only model optically thin
emission. The linear spectral decomposition gives erroneous results if in some
spectral regions the emission is optically thick. The latter constrain
motivated our choice to restrict the fitting with `specfit` to the $N$ band.
To overcome that limitation, we apply our own model fitting tool,
TGMdust202020TGMdust is available for download at
https://home.strw.leidenuniv.nl/~varga/pro/tgmdust. (Temperature Gradient
Modeling with dust), which is a python implementation of our multi-zone disk
model described in Sect. 3. For the optimization we use `emcee`, an MCMC
ensemble sampler (Foreman-Mackey et al. 2013a, b). In `TGMdust` we use our
original assumption for the radiative transfer, described in Eq. 4, capable to
represent optically thick emission. We note that the linear combination
assumption employed in the previous step provides nearly equivalent results to
the `TGMdust` fits if the dust emission is optically thin. In this step of
modeling, we use all data we have, not just the $N$ band. The large number of
model parameters (44 in the case of the 3-zone model) might prevent the MCMC
sampler from reaching convergence. Thus, with TGMdust we fit a limited set of
parameters, and we use the dust mass fractions obtained in the first step by
specfit as fixed parameters.

##  Appendix F Degeneracy in the spectral decomposition

As was shown by van Boekel et al. (2005) (see Sect. 5.2 and Fig. 8 therein),
the silicate spectral templates used in the spectral decomposition are
degenerate to a degree. We investigated this issue in our specfit model run by
inspecting the individual samples. We found several good-fitting solutions
with a dust composition differing considerably from the best fit. One example
is shown in Table 6, and in Figures 18, 19. The $\chi^{2}$ of the alternative
solution is $1.7$ times larger than that of the best fit. The fit is
qualitatively similar to that of the best fit, except at the edges of the band
($<8.5\ \mu$m and $>12.5\ \mu$m) where there is a larger deviation from the
data. The zone boundaries ($0.2$, $0.9$, and $4.0$ au) are close to the values
of the best fit, while the mass fraction values indicate degeneracies between
large pyroxene and enstatite in zone 2, and between large olivine, small
pyroxene, and enstatite in zone 3. These findings imply that despite the
degeneracies in the silicate spectral templates, the disk substructrures (zone
boundaries) are robustly constrained by our modeling.

Table 6:  An example for an alternative good-fitting solution for the silicate spectral decomposition. The abundances are given as fractions of the total dust mass in percent, excluding the hidden dust component responsible for the continuum emission. Species | Zone 1 | Zone 2 | Zone 3
---|---|---|---
|  $0.2-0.9$ au |  $0.9-4.0$ au |  $>4.0$ au
$\text{Olivine}\ 2.0\ \mu\text{m}$ | $14.9$ | $38.7$ | $0.0$
$\text{Pyroxene}\ 0.1\ \mu\text{m}$ | $0.0$ | $43.4$ | $50.9$
$\text{Pyroxene}\ 2.0\ \mu\text{m}$ | $9.3$ | $15.8$ | $0.0$
$\text{Forsterite}\ 0.1\ \mu\text{m}$ | $2.9$ | $1.7$ | $1.2$
$\text{Enstatite}\ 5.0\ \mu\text{m}$ | $72.9$ | $0.4$ | $47.9$
Crystallinity fraction | $75.8$ | $2.1$ | $49.1$

Figure 18: Same as Fig. 1, but for a specfit solution corresponding to the
dust mass fractions listed in Table 6. Figure 19:  Same as Fig. 2, but for a
specfit solution corresponding to the dust mass fractions listed in Table 6.

##  Appendix G Supplementary model plots

###  G.1 TGMdust model run with iron grains

Figure 20: Corner plot showing the posterior distributions of the MCMC samples
from the run with iron grains. We note that in case of the surface densities
and mass fractions we fit the logarithm of the quantities.

###  G.2 TGMdust model run with carbon grains

Figure 21: Best-fit model images at various data wavelengths with the TGMdust
run with carbon grains. The brightness scaling is linear and homogeneous
across all wavelengths. For a better perception of the fainter structural
features, the images become saturated at $3.4\times 10^{-12}$ W m-2 Hz-1 sr-1.
The gray circles in the bottom left corner indicate the approximate beam size
(estimated as $0.77\lambda/(B_{\text{p,max}})$, where $B_{\text{p,max}}$ is
the maximum baseline length). Figure 22: Radial surface brightness profiles of
the disk corresponding to the best-fit TGMdust model with carbon grains. The
shaded area around each curve is the $16-84$ percentile range from the last
$5000$ samples of the chain.

Figure 23: Vertical optical depth of the best-fit disk model with carbon
grains. Left panel: Radial optical depth profiles at selected wavelengths.
Right panel: Optical depth spectra at the inner edges of each zone.

Figure 24:  Corner plot showing the posterior distributions of the MCMC
samples from the run with carbon grains. We note that in case of the surface
densities and mass fractions we fit the logarithm of the quantities.