\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,published]{cedram}
\newenvironment{noXML}{}{}%
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}}
\def\tabnote#1{\vskip4pt\parbox{.935\linewidth}{#1}}
\usepackage{etoolbox}
\usepackage{upgreek}
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\def\xsection#1{}
\newenvironment{inftab}{}{}
\newcounter{runlevel}
\usepackage{multirow} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20240924}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\0{\phantom{0}}
\def\ndash{\text{--}}
\def\rmmu{\upmu}
\def\mn{\phantom{$-$}}
\def\xsmath#1{\ensuremath{#1}}
\def\refinput#1{}
\def\back#1{}
\def\r#1{\mathring{#1}}
\def\xensuremath#1{\ensuremath{#1}}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\DOI{10.5802/crmeca.275}
\datereceived{2024-11-07}
\dateaccepted{2024-11-19}
\ItHasTeXPublished
\newtheorem{theorem}{Theorem}%
\renewcommand\thetheorem{\unskip}
}
{\documentclass[crmeca]{article}
%\usepackage[T1]{fontenc}
\usepackage{upgreek}
\usepackage{stmaryrd}
\def\CDRdoi{10.5802/crmeca.275}
\def\r#1{\mathring{#1}}
\def\xsmath#1{${#1}$}
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
}
\makeatother

\usepackage{upgreek}

\dateposted{2025-01-08}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{review}

\title{Gluing for the Einstein constraint equations}

\alttitle{Recollement pour les \'equations de contrainte d'Einstein}

\author{\firstname{Justin} \lastname{Corvino}\CDRorcid{0009-0001-6914-1253}}
\address{Lafayette College, Department of Mathematics, Easton, PA 18042, USA}
\email{corvinoj@lafayette.edu}

\keywords{\kwd{Einstein constraint equations}\kwd{Geometric
analysis}\kwd{Gluing constructions}}

\altkeywords{\kwd{Equations de contrainte d' Einstein}\kwd{Analyse
g\'eom\'etrique}\kwd{M\'ethodes de recollement}}

\thanks{AMS-Simons Research Enhancement Grant for PUI Faculty}

\begin{abstract} 
Initial data for the Einstein equations must satisfy a system of
nonlinear partial differential equations, the Einstein constraint
equations. Constructing interesting solutions of the constraint
equations can then lead to interesting spacetime evolutions. Over the
past twenty-five years, various gluing methods to construct solutions
of the constraints have shed light on some long-standing questions. In
this article we survey some of the methods and results achieved to
date.
\end{abstract} 

\begin{altabstract} 
Les donn\'ees initiales des \'equations d'Einstein doivent satisfaire
un syst\`eme d'\'equations aux d\'eriv\'ees partielles non lin\'eaires,
les \'equations de contrainte d'Einstein. La construction de solutions
int\'eressantes des \'equations de contrainte peut alors conduire \`a
des \'evolutions int\'eressantes de l'espace-temps. Au cours des
vingt-cinq derni\`eres ann\'ees, diverses m\'ethodes de recollement
permettant de construire des solutions des \'equations de contrainte
ont permis de faire la lumi\`ere sur certaines questions de longue
date. Dans cet article, nous passons en revue quelques-unes des
m\'ethodes et des r\'esultats obtenus \`a ce jour.
\end{altabstract}

\maketitle
{\vspace*{-.6pc}}

\end{noXML}

\section{Introduction and preliminaries}

The foundational work of Yvonne Choquet-Bruhat\footnote{For her
exposition on this, see her book~\cite{ycb:grein}, the classic survey
articles~\cite{cb-cp, cb-yo}, as well as~\cite{CB-BCP} for her comments
on some of the history of the problem.} firmly established that the
Einstein field equation can be treated as an evolution
problem~\cite{CB52}. The initial data for the vacuum equation are given
by a Riemannian manifold $(M, g)$ along with a symmetric $(0,2)$-tensor
$K$, that determine a vacuum spacetime, i.e.\ a time-orientable
Ricci-flat Lorentzian manifold $(\mathcal{S}, \bar{g})$ in which $M$ is
a Cauchy surface embedded with first and second fundamental forms $g$
and $K$, which play the role of the initial metric and its time
derivative (in appropriate coordinates). In the non-vacuum case, the
initial data will be augmented with data for whatever physical fields
are being modeled. 

The initial data $g$ and $K$ are not freely prescribed, but must adhere
to the Gauss and Codazzi equations relating the ambient geometry to the
fundamental forms. The Einstein equation involves the Ricci tensor, and
so from appropriate traces of the Gauss and Codazzi equations along
with the Einstein equation, one deduces a coupled, nonlinear {\it
underdetermined} elliptic system for $g$ and $K$, the exploration of
which has provided a rich playground for mathematical physicists and
geometric analysts, with deep connections to topology and geometry, in
part forged by advances in the theory of partial differential
equations. Robust methods to construct solutions of the constraint
equations allow for constructions, theoretical and numerical, of
spacetimes. Understanding the solution space to the constraint
equations is thus intimately tied to an understanding of the solution
space of the Einstein field equation, and a parametrization of the
solution space of the constraints would ostensibly isolate the true
gravitational degrees of freedom~\cite{yo:gf}, cf.~\cite{max-drift}. 

A fulsome account of the work on the constraint equations, which one
can trace back at least eighty years to classical work of
Lichnerowicz~\cite{lich-ce}, is well beyond the scope of this article. 
Rather, we focus on some relatively recent work of the past twenty-five
years on constructing solutions of the Einstein constraint equations by
{\it gluing} techniques, discussing some of the methods and
applications without giving an exhaustive account. Gluing techniques
can be used to address the following natural question: given multiple 
solutions of the constraint equations, is there another solution (on a
connected manifold) which contains chosen regions from the given
initial data sets (at least approximately)? There is of course the
related question about spacetimes: given regions of multiple spacetimes
solving Einstein's field equation, is there a (connected) spacetime
satisfying the field equation that contains them? 

From the initial data point of view, if one can construct an initial
data set which contains regions obtained from several spacetimes, the
resulting evolution of the initial data will contain regions isometric
to ones in the given spacetimes, by finite speed of propagation. That
said, one might naturally try to glue spacetime regions together
directly, or at least utilize the Lorentzian structure in the
construction. In this direction we mention two different recent
series of works. Hintz constructs spacetimes with multiple small black
holes by a gluing construction using geometric singular analysis
(pseudodifferential $b$-calculus) techniques (see~\cite{Hin23},
cf.~\cite{Hin22}). In a different direction, Aretakis, Czimek and
Rodnianski achieve a suite of gluing results by {\it characteristic
gluing}, which leverages the {\it null structure} of spacetime, as
surveyed with additional references
in~\cite{aretakis-czimek-rodnianski}. While these works are
interesting, we do not have the time in the present short survey
article to do justice to either, and we will focus on constructing
solutions of the constraint equations on a spacelike hypersurface
(though we do mention some compelling results obtained by the methods
of~\cite{aretakis-czimek-rodnianski} in Section~\ref{sec:further}
below). 

\subsection{Overview of some spacelike gluing results} \label{sec1.1}

A natural interpretation of the above aim of gluing is that of
constructing initial data for a relativistic $N$-body
problem~\cite{cci-N-cqg}, including, say, multi-black hole spacetimes. 
Such constructions have been achieved in earlier works~\cite{ptc-maz,
cci-ts, cci}, while more recent works have used gluing to construct
more exotic such initial data sets, such as the $N$-body configurations
with gravitational shielding from~\cite{car-sch}, multi-localized
initial data sets~\cite{acp:multi}, as well as initial data with small
black holes~\cite{Hin22}. Naturally, the issue of the two-body and
$N$-body problems in general relativity dates back a century
(see~\cite{rac-th, darm, lich-ce} and references therein) and various
sorts of approaches have been taken to produce initial data to model
such spacetimes. 

\looseness=-1
Connected sum construction of initial data was achieved by Isenberg,
Mazzeo and Pollack in~\cite{imp} by the conformal method, and
Chru\'sciel and Delay~\cite{cd} observed how to combine conformal and
localized gluing techniques (introduced in~\cite{cor:schw, cs:ak}) to
perform connected sum constructions which leave the data unperturbed
outside a connecting neck region (cf.~\cite{cip, cem} for related
results, and~\cite{mazz:surg-imp, del-mazz:ref-glue} for analogous
surgery results). Of particular note,~\cite{imp:top} established that
there is no topological obstruction to admitting an asymptotically
Euclidean solution of the vacuum Einstein constraint equations: for a
closed manifold $\Sigma$ and $p\in \Sigma$, there is a vacuum
asymptotically flat solution of the Einstein constraints on $\Sigma
\setminus \{p\}$. An intriguing implementation of the conformal
gluing method was employed by Stavrov in~\cite{stav-ppl} to construct
families of initial data to model point particle limits in general
relativity (cf.\ recent work of Hintz~\cite{Hin22, Hin23} along these
lines). 

Gluing constructions have also shed light on the asymptotics of
solutions of the constraints. The original motivation for the
construction in~\cite{cor:schw} was to understand to what extent the
asymptotics determine the interior of a vacuum initial data set, in
particular how special are Schwarzchild asymptotics. That initial
work showed that any compact subset of any time-symmetric
asymptotically flat vacuum initial data set is also contained in a
vacuum time-symmetric initial data set which is {\it exactly}
Schwarzschild near infinity. While that in itself may seem intriguing,
one immediately turns to the spacetime evolution of such data. To be
precise, the question open at the time was whether there exist any
nontrivial asymptotically simple vacuum spacetimes (sometimes called
{\it purely radiative} spacetimes), a question answered in the
affirmative in~\cite{cd:as-pen, cor:as-pen} by constructing suitable
families of initial data which are Schwarzschild with small mass in a
{\it uniform} neighborhood of infinity, evolving to an asymptotically
hyperboloidal hypersurface, and then invoking the work
of~\cite{friedrich} for the asymptotics of the evolution. 

Another question that has been addressed by gluing methods is that of
{\it gravitational shielding}. A remarkable result was obtained by
Carlotto and Schoen~\cite{car-sch}, who showed that any region of any
vacuum initial data set can be localized within a conical region in a
vacuum initial data set which is {\it flat} outside the conical region.
While Bartnik's parabolic quasi-spherical construction~\cite{ba:qs}
produced examples of nontrivial asymptotically flat vacuum initial data
(time-symmetric, so scalar-flat) containing a flat region, their
localized gluing construction produced striking examples with
non-compact flat regions. The results in~\cite{car-sch} were proven in
spirit along the lines of~\cite{cor:schw, cs:ak, cd}, with a number
additional of technical obstacles to surmount. While the
works~\cite{cor:schw, cs:ak, cd} established that any bounded region in
an asymptotially flat vacuum solution can be realized as a subdomain in
a solution in which each asymptotic end is 
a slice of a Schwarzschild or Kerr spacetime, the method of proof
requires that the end be tuned appropriately to the given data (in
terms of the asymptotically conserved quantities: ADM energy-momentum,
center-of-mass and angular momentum), whereas the regions in the result
of Carlotto and Schoen are glued to flat exteriors outside the conical
regions. In any case, all these results illustrate that unique
continuation fails spectacularly for the vacuum constraints.

In this article we overview aspects of gluing by the conformal
method~\cite{imp} and compare with the localized approach. While there
are many excellent references for the conformal method, localized
gluing constructions are relatively newer, and are still being
developed, so we will focus on some details of the approach to
localized gluing that dates back to the works~\cite{cor:schw, cs:ak}. 
In the spirit of ``something old and something new'', we illustrate the
localized gluing technique by applying it to the time-symmetric
Einstein--Maxwell constraint equations.

{\vspace*{-3pt}} 
\subsection{Preliminaries} \label{sec1.2}
\subsubsection{The Einstein constraint equations} \label{sec:ece}

The vacuum constraint equations arise from the following. Consider a
spacetime splitting $(\mathcal{S} =I\times M, \bar{g})$, with $I$ an
interval and $(M, g)$ a Riemannian hypersurface 
given as a level set of the coordinate $x^0$
on $I$ (so ${\partial}/{\partial x^0}$ is timelike). With $G=
\mathrm{Ric}(\bar{g}) - (1/2) R(\bar{g}) \bar{g}$ the Einstein tensor,
$G_{\alpha \beta} n^\beta $ does not contain any second time
derivatives $\bar{g}_{\mu \nu, 00}$ of the metric components, where 
$n$ is a timelike unit normal to $M$. In fact, using the Gauss and
Codazzi equations, one obtains, with $G(n, \cdot)$ a one-form along $M$, 
{\begin{eqnarray}
& G(n, n) = \tfrac12 (R(g) - |K|^2_g + (\mathrm{tr}_g K)^2) \label{Gnn}
\Seqnsplit
& G(n, \cdot) = \mathrm{div}_g (K- (\mathrm{tr}_g K)g). \label{Gn0}
\end{eqnarray}}\unskip
If we consider a cosmological constant $\Lambda$, and let $G_{\Lambda}
= G+ \Lambda \bar{g}$, and if $T$ is the stress-energy tensor, the
Einstein equation takes the form $G_{\Lambda} = \kappa T$ for some
constant $\kappa>0$. (In spacetime dimension four, $\kappa = {8\uppi
G}/{c^4}$.) We observe that $\mu := T(n,n)$ is the energy density as
measured by an observer with four-velocity $cn$, and $c^{-1} J :=
-c^{-1}T(n, \cdot)$ is corresponding momentum density one-form. The
dominant energy condition that the energy-momentum vector of the matter
fields is causal future-pointing is given by $\mu\geq |J|_g$. Combining
the Einstein equation with (\ref{Gnn})--(\ref{Gn0}) yields the Einstein
constraint equations
{\begin{eqnarray}
& \tfrac12(R(g) - |K|^2_g + (\mathrm{tr}_g K)^2) = \kappa \mu +
\Lambda \label{ham}\Seqnsplit
& \mathrm{div}_g (K- (\mathrm{tr}_g K)g) = \kappa J. \label{mom} 
\end{eqnarray}}\unskip

In the vacuum ($T=0$) time-symmetric ($K=0$) case, the system reduces
to the constant scalar curvature (CSC) condition $R(g) = 2\Lambda$. 
Unless noted, we take $\Lambda =0$, in which case the time-symmetric
vacuum constraints reduce to $R(g)=0$.

For a non-vacuum example, consider initial data for the time-symmetric
Einstein--Maxwell equations, often given as a Riemannian metric $g$ and
an electric field~$E$. The constraint (\ref{ham}) becomes $R(g) =
2\kappa |E|^2_g$. In any source-free region, we must also have
$\mathrm{div}_g E= 0$. For simplicity we will treat $E$ as a vector
field and study the map $\Psi(g, E) = (R(g) - 2\kappa |E|_g^2,
\mathrm{div}_g E)$. We note that it could be more natural to specify
the field through a differential form $\eta_E$, to which a metric $g$
would associate a vector field $E_g$ in such as way that the condition
$\mathrm{div}_g E_g=0$ would follow {\it for any} $g$ in a source-free
region. We will see the echo of this in the gluing method below. For
further exploration of specifying the initial data for fields,
see~\cite{isen-max-phase}. 

We will focus on the constraints in dimension $n=3$, though much of
what is discussed below extends readily to $n\geq 3$. We employ the
Einstein summation convention of summing over repeated upper and lower
indices; a semicolon denotes a covariant derivative, and a comma
denotes a partial derivative. 

\subsubsection{Some classical solutions} 

For later use, we recall a family of time-symmetric solutions to the
Einstein--Maxwell constraint equations, each obtained from data on a
spacelike hypersurface in a Reissner--Nordstr\"{o}m spacetime: 
$g_{\mathrm{RN}}= (1- ({2m}/{r}) + ({Q^2}/{r^2}))^{-1}\, \mathrm{d}r^2+
r^2 \r{g}_{\mathbb{S}^2}$, with $E_{\mathrm{RN}}= ({Q}/{r^2})
\sqrt{1- ({2m}/{r}) + ({Q^2}/{r^2})} ({\partial}/{\partial r}) =:
({Q}/{r^2}) e_r$, where $\r{g}_{\mathbb{S}^2}$ is the unit round
metric. Since $R(g_{\mathrm{RN}}) = {2 Q^2}/{r^4}$, we see the units
are such that $\kappa=2$. We will just employ the data on an
asymptotic end $r=|x|\gg 1$, and so we do not pose any relation between
$m$ and $Q$, though that could affect the Reissner--Nordstr\"{o}m
spacetime structure.

For $Q=0$, the Reissner--Nordstr\"om metric reduces to the
Schwarzschild metric, which has vanishing scalar curvature. We write
the metric in isotropic form, and we can also re-center the coordinates
to obtain $(g_S^{m, c})_{ij}(x)=(1+ {m}/({2|x-c|}))^4 \delta_{ij}$. A
metric $g_{ij}=u^4 \delta_{ij}$ which is conformally Euclidean has
vanishing scalar curvature precisely when $u$ is Euclidean harmonic. 
If we require that $u$ tends to $1$ at infinity, to impose
asymptotically Euclidean geometry, then $u$ is either identically 1 or
it must have singularities, such as for the Schwarzschild metric with
$m\neq 0$, and for the following family of solutions. 

As recalled earlier, it has been a longstanding question as to how to
construct solutions to the Einstein equations to model the interaction
of multiple bodies. A classical solution of the vacuum constraints is
given by vacuum Brill--Lindquist metrics (see~\cite{bl, Sor-Stav},
cf.~\cite{lich-ce}): given a finite set $\mathfrak{p}=\{p_1,\ldots,
p_N\}$ of $N$ points in $\mathbb{R}^3$ and $\mathfrak{m}=(m_1, \ldots,
m_N)\in \mathbb{R}^N$, we consider the metric
{\begin{equation} \label{eq:BLm}
(g^{\mathfrak{m}, \mathfrak{p}}_{\mathrm{BL}})_{ij} (x) =
\left(1+ \sum_{k=1}^N \frac{m_k }{2|x-p_k|}\right)^4 \delta_{ij}
\end{equation}}\unskip
for $x\in \mathbb{R}^3 \setminus \mathfrak{p}$. This metric is
conformally Euclidean with harmonic conformal factor, and thus has
vanishing scalar curvature, providing a time-symmetric vacuum initial
data set. If we take all $m_k>0$, then the data has $N+1$
asymptotically flat ends. By an ingenious use of the method of images,
Misner~\cite{Mis63} modified the conformal factor, incorporating an
infinite sequence of inversions, to produce a scalar-flat metric with
two asymptotically flat ends and $N$ Einstein--Rosen bridges, hence
obtaining a type of gluing construction by conformal means. In another
direction, one can modify the Brill--Lindquist data by localized gluing
techniques to produce data on $\mathbb{R}^3$, even allowing
$\mathfrak{p}$ to be infinite; see~\cite{acp:multi} along with comments
in Section~\ref{sec:further} below. 

\subsubsection{The conformal method} \label{sec1.2.3}

The vacuum Einstein constraints form an {\it underdetermined} elliptic
system for $g$ and $K$. The conformal method specifies some parts of
$g$ and $K$, leaving a determined elliptic system. This approach can
be traced back to\footnote{One can find an early analysis of the
conformally Euclidean case in~\cite{rac-th}, which also includes some
early study of static spaces and asymptotically Euclidean spaces; 
see~\cite{CB-BCP} for more comments on the history.}
Lichnerowicz~\cite{lich-ce}, Choquet-Bruhat~\cite{cb-jrm}, York and \'O
Murchadha~\cite{yo:gf, om-yo-73, om-yo-prd1}. The formulation below is
especially well suited to the constant mean curvature (CMC) case with
$\mathrm{tr}_g K$ constant\footnote{See~\cite{max-CTS},
cf.~\cite{max-drift, isen-max-phase}, for discussion of the conformal
method outside the vacuum CMC regime.} (for a variation
see~\cite{om-yo-prd1}). It is based on the transverse-traceless (TT)
decomposition of symmetric two-tensors. In this approach, a conformal
class $\mathcal{C}$ of metrics is given, say $\mathcal{C}= [\gamma]$. 
Furthermore a scalar function $\tau$ and a $(0,2)$-tensor $\sigma$ are
given on $M^n$, where $\sigma$ is TT with respect to $\gamma$:
$\mathrm{div}_{\gamma} \sigma =0$, $\mathrm{tr}_{\gamma} \sigma =0$. 

From this data we seek a scalar function $u>0$ and a vector field $W$
such that, with
{$$
g= u^{\frac{4}{n-2}} \gamma, \quad K= u^{-2} (\sigma + L_{\gamma} W)
+ \frac{\tau}{n} u^{\frac{4}{n-2}} \gamma,
$$}\unskip
where $(L_{\gamma} W)_{ij} = W_{i;j} + W_{j;i} - ({2}/{n})
(\mathrm{div}_{\gamma} W)\gamma_{ij}$ is the conformal Killing
operator, $(g, K)$ satisfies the vacuum constraint equations, which
take the form
{\begin{eqnarray}
&\displaystyle \Delta_{\gamma} u - \frac{n-2}{4(n-1)} R(\gamma) u +
\frac{n-2}{4(n-1)} |\sigma + L_{\gamma} W|^2_{\gamma} u^{-\frac{3n-2}
{n-2}} - \frac{n-2}{4n} \tau^2 u^{\frac{n+2}{n-2}}=0 \label{eq:hccm}
\Seqnsplit
&\displaystyle \mathrm{div}_{\gamma} (L_{\gamma} W) - \frac{n-1}{n}
u^{\frac{2n}{n-2}}\,\mathrm{d}\tau =0. \label{eq:mccm}
\end{eqnarray}}\unskip
The system is a semilinear elliptic system for $(u, W)$. In the CMC
case (constant~$\tau$), the system de-couples: the second equation
becomes $\mathrm{div}_{\gamma} (L_{\gamma} W) = 0$. On a closed
manifold, or asymptotically flat manifold with $W$ decaying suitably
near infinity, multiplying by $W$ and integrating by parts gives
$L_{\gamma} W=0$, i.e.\ $W$ is a {\it conformal Killing field}. In this
case $W$ does not affect $K$ or the equation for~$u$, and we can just
take $W=0$ (which is the only option in the absence of conformal
symmetries). Then the conformal method comes down to solving the {\it
Lichnerowicz equation}
{$$
\Delta_{\gamma} u - \frac{n-2}{4(n-1)} R(\gamma) u + \frac{n-2}
{4(n-1)} |\sigma|^2_{\gamma} u^{-\frac{3n-2}{n-2}} - \frac{n-2}{4n}
\tau^2 u^{\frac{n+2}{n-2}} = 0.
$$}\unskip
As noted earlier, we will focus on $n=3$, in which case the above
becomes
{\begin{equation} \label{eq:lichn3}
\Delta_{\gamma} u - \tfrac18 R(\gamma) u + \tfrac18
|\sigma|^2_{\gamma} u^{-7} - \tfrac{1}{12}\tau^2 u^5=0.
\end{equation}}\unskip

\subsubsection{The conformal method and asymptotics} \label{sec1.2.4}

When gluing initial data with the conformal method (see
Section~\ref{sec:con-glue}), the resulting data can be made arbitrarily
close to the given data outside the gluing region. A similar scenario
arises when using conformal perturbations to arrange conformally
Euclidean, or more generally {\it harmonic}, asymptotics. Indeed Schoen
and Yau~\cite{sy-elm} showed that a scalar-flat asymptotically flat
metric can be approximated by one which is conformally Euclidean
outside a compact set, and so that the ADM energy is perturbed by
an arbitrarily small amount. The method of proof is straightforward:
along a large annular region, patch the original metric to a Euclidean
metric near infinity; for large enough annuli, we get approximately
scalar-flat data. One then employs a conformal factor to reimpose the
vacuum constraint, with appropriate estimates on the size of the
perturbation of the metric and the ADM energy. Such an approximation
result is useful for proving the Positive Energy Theorem~\cite{sy:pmt}
and the Riemannian Penrose Inequality, since it simplifies the
asymptotics with only a small change to the energy, and Bray termed
such metrics {\it harmonically flat} at infinity~\cite{bray-pen-jdg}. 

An analogous non-time-symmetric version established in~\cite{cs:ak},
{\it harmonic asymptotics}, has also proved useful in this and other
contexts~\cite{ehls-pmt, huang-cm-cqg}. The formulation is a
modification of the conformal method above so that the principal part
of the constraints operator in the asymptotic region is a diagonal
Laplacian. The asymptotic form of the data in harmonic asymptotics is
{\begin{equation} \label{eq:harm-asy}
g_{ij} = u^4 \delta_{ij}, \quad K_{ij} = u^2 \left(X_{i, j} + X_{j,i}
- \tfrac12 X^k_{,\,k}\delta_{ij}\right).
\end{equation}}\unskip
This form is usually written in terms of the momentum tensor $\uppi=K -
(\mathrm{tr}_g K) g$, in which case this becomes $\uppi_{ij} = u^2
(X_{i, j} + X_{j,i} - X^k_{,\,k}\delta_{ij})$.

The proof for harmonic asymptotic approximation of asymptotically flat
initial data sets is similar to the time-symmetric case. Given
asymptotically flat data $(g, K)$, smoothly patch the metric to
Euclidean while cutting off $K$ to zero, over a large annular region
$R<|x|<2R$, with the resulting data $(\hat{g}, \hat{\uppi})$ (a family
of data depending on~$R$). Let $(\mathcal{L}_g X)_{ij} =
X_{i;j}+X_{j;i} - (\mathrm{div}_g X)g_{ij}$. One seeks to reimpose the
constraints and obtain harmonic asymptotics by considering data of the
form $\bar{g}= u^4 \hat{g}$ and $\bar{\uppi} = u^2 (\hat{\uppi} +
\mathcal{L}_{\hat{g}} X)$, or $\bar{K} = u^2 (\hat{K}
+\mathcal{L}_{\hat{g}} X + (1/2) (\mathrm{div}_{\hat{g}} X)\hat{g})$.
Observe that it is not the conformal Killing operator $L_{\hat{g}} X$
that appears here, in contrast to the standard formulation of the
conformal method. In any case, $(\hat{g}, \hat{K})$ is an approximate
solution of the vacuum constraints (better for larger~$R$), and one
seeks to perturb to an exact solution, for which one considers the map
{$$
(u, X) \mapsto T(u, X)=\left(R(\bar{g}) - |\bar{K}|^2_{\bar{g}} +
(\mathrm{tr}_{\bar{g}} \bar{K})^2, \mathrm{div}_{\bar{g}} (\bar{K}-
(\mathrm{tr}_{\bar{g}} K)\bar{g})\right).
$$}\unskip
Outside a compact set, i.e.\ for $|x|>2R$, the vacuum constraints $T(u,
X)=0$ can be written in Cartesian components as (omitting metric
subscripts for the Euclidean metric)
{\begin{equation} \label{eq:harm-asy-con}
8 \Delta u = u \left({-}|\mathcal{L} X|^2 + \tfrac12 (\mathrm{tr}
(\mathcal{L} X))^2\right),\quad \Delta X^i=-4u^{-1}u_{,j}(\mathcal{L}
X)^j_{i} +2 u^{-1} u_{, i} \mathrm{tr}(\mathcal{L} X).
\end{equation}}\unskip
Solutions in spaces with suitable decay for $u-1$ and $X$ will admit
partial expansions of the form $u(x) = 1+ ({a}/{|x|}) + O(|x|^{-2})$,
and $X^i(x) = ({b^i}/{|x|}) + O(|x|^{-2})$. Using (\ref{eq:harm-asy}),
one finds the ADM energy is $2a$, and the linear momentum is $P^i =
-(1/2) b^i$.

The linearization of $T(u, X)$ at $(1,0)$ is readily found to be 
Fredholm of index zero between appropriate spaces~\cite{cs:ak}, of
which (\ref{eq:harm-asy-con}) is certainly indicative. The possible
existence of a finite-dimensional cokernel is accommodated by finding a
complementary finite-dimensional space of {\it compactly supported}
tensors, and showing that one can find $(h, k)$ near $(0,0)$ in this
complementing space along with $(u, X)$ close to $(1,0)$ (for large
$R$), such that $(\bar{g} + h, \bar{K} + k)$ solves the vacuum
constraints. It can also be shown that the energy and linear momentum
can be perturbed an arbitrarily small amount in this
construction~\cite{cs:ak}. 

In a certain sense, these approximation constructions establishing that
a kind of asymptotic behavior is suitably dense can be construed as a
type of gluing construction, gluing a given interior to an exterior end
in a certain family. Good approximate solutions form a family
parametrized by~$R$, which can be corrected with suitable
perturbations. In this case the family of exterior data, either the
harmonically flat metrics, or more generally data with harmonic
asymptotics, is a large family, as compared to a data like the
Schwarzschild family, parametrized by $m$ and~$c$. Just as in the
conformal gluing we recall below, the given data is perturbed a small
(for large~$R$) amount. 

Harmonic asymptotics yield a particularly convenient form of the
asymptotics for the energy and momenta. For some other applications,
one might want to make controlled construction of solutions of the
constraints with more general asymptotics. For instance, recent work of
Fang, Szeftel and Touati~\cite{FST-mink, FST-BH} constructs suitable
classes of vacuum initial data compatible with the stability results
for Minkowski spacetime and for black hole stability. The basic
approach is to start with an approximate solution with a desired
asymptotic behavior and then perturb to a solution of the constraints,
with estimates on the perturbation which ensure the desired decay rate.
The perturbation is as above, utilizing a conformal change of metric,
together with a perturbation of the momentum tensor of the form
$\mathcal{L}_{\gamma} X$ (for a suitable metric~$\gamma$), plus a
complementing compactly supported perturbation. 

\subsubsection{Remarks on localized deformations} \label{sec1.2.5}

A fundamental issue for constructing solutions to the constraints by
gluing methods is how to obtain controlled perturbations of the
constraints map (where $\uppi = K - (\mathrm{tr}_g K) g$)
{\begin{equation} \label{eq:vac-Phi}
(g, \uppi) \mapsto \Phi(g, \uppi) := (R(g) - |K|^2_g +
(\mathrm{tr}_gK)^2, \mathrm{div}_g \uppi).
\end{equation}}\unskip
As recalled above, modifying data with a suitable conformal
transformation has proven fruitful. Since a conformal factor $u$
solves an elliptic equation, we do not expect to be able to effectively
contain the support of $u-1$. With the linearization stability results
of Fischer and Marsden~\cite{fm-sc-def, fm:defcon} in mind, one can
naturally move outside the conformal regime and consider a general form
of perturbations, modifying $(g, \uppi)$ to $(g+\delta g, \uppi+\delta
\uppi)$. The vacuum constraints then give an underdetermined elliptic
system for the perturbations $(\delta g, \delta \uppi$). This freedom
is what has been employed in~\cite{cor:schw, cs:ak} and numerous
subsequent works to achieve more direct control on the support of the
deformation tensors. 

Of course, the analysis will focus on the linearized constraints
operator, which morally speaking should be locally surjective when the
formal adjoint of the linearization has trivial kernel. The presence
of nontrivial kernel indicates symmetries in the spacetime evolution of
the data, as in~\cite{monc:stat}, so that except in special situations,
we expect to be able to effectively perturb the constraints operator. 
That said, in many gluing constructions, the approximate solution
approaches a special configuration which admits symmetries, and this
requires some extra care to handle. 

To make this precise, we consider the linearization and formal
$L^2$-adjoint for the scalar curvature and vacuum constraints operators.
The well-known formulas for the linearization of the scalar curvature
operator and its adjoint are given by 
{\begin{equation}\label{eq:scal-lin}
DR_g(h) = -\Delta_g (\mathrm{tr}_g h)+ \mathrm{div}_g \mathrm{div}_g
h - h \cdot_g \mathrm{Ric}(g), \quad DR_g^*(f) = - (\Delta_g f) g +
\mathrm{Hess}_g f - f\mathrm{Ric}(g).
\end{equation}}\unskip
At a Ricci-flat metric, an element of the kernel of $DR_g^*$ is
precisely a function $f$ with vanishing Hessian, which in the Euclidean
case means $f$ is a constant plus a linear combination of the Cartesian
coordinate functions $x^1, x^2, x^3$. 

While the formulas for the linearization of $\Phi$ and its formal
adjoint are straightforward to compute in the general case, they
simplify tremendously at the flat data, about which we will be
perturbing. Let \xsmath{\r{g}} be the Euclidean metric
on~$\mathbb{R}^3$. If we let \xsmath{D\r{\Phi} = D\Phi_{(\r{g}, 0)}},
we find from (\ref{eq:vac-Phi})
{$$
D\r{\Phi}(h, \omega) = (DR_{\r{g}} (h),
\mathrm{div}_{\r{g}} \omega), \quad D\r{\Phi}^*
(f, Z) = \left(DR^*_{\r{g}} (f), -\tfrac12 L_Z \r{g}\right),
$$}\unskip
where $(L_Z g)_{ij} = Z_{i;j} + Z_{j;i}$ is the Lie derivative, so that
at \xsmath{\r{g}} in Cartesian coordinates, we have \xsmath{(L_Z
\r{g})_{ij} = Z_{i,j} + Z_{j,i}}. Now, $L_Z g=0$ precisely for a
Killing vector field~$Z$, which at \xsmath{\r{g}} must be a linear
combination of the generators of translations and rotations. Thus the
kernel of \xsmath{D\r{\Phi}^*} gives a ten-dimensional space of
potential obstructions to localized gluing constructions for the vacuum
Einstein equations, near Minkowskian initial data. To accommodate for
this, we glue to model families of solutions that are governed by ten
parameters which can cover this approximate cokernel, such as the Kerr
family of asymptotically flat initial data (see, e.g.\ 
\cite[Appendix~F]{cd}). We will see in Section~\ref{sec:loc-glue} below
how the mass and center-of-mass parameters can be used to handle the
kernel of \xsmath{DR_{\r{g}}^*}, while in the general case, the linear
and angular momentum parameters correspond to the translational and
rotational Killing fields~\cite{cs:ak, cd}. In general the kernel
elements of $D\Phi^*_{(g, \uppi)}$ are related to spacetime
symmetries~\cite{monc:stat}, and are called {\it Killing Initial
Data}.\looseness=-1

\section{Gluing by the conformal method} \label{sec:con-glue}

We now discuss the connected sum construction of Isenberg, Mazzeo and
Pollack~\cite{imp}. We start with solutions of the vacuum constraints
on manifolds $M_1$ and $M_2$ with constant mean curvature~$\tau$, given
on each as a Riemannian metric~$\gamma_j$, and a TT tensor $\mu_j$ with
respect to~$\gamma_j$, so that the second fundamental form is given
$K_j= \mu_j + ({\tau}/{3}) \gamma_j$ (recall that we consider dimension
three). 

To construct a connected sum $M=M_1\sharp M_2$ of disjoint manifolds
$M_1$ and $M_2$, one removes a ball around points $p_1\in M_1$ and
$p_2\in M_2$ and identifies the union of the complements along the
boundary spheres (the construction could also be applied where $p_1$
and $p_2$ are points on the same manifold, corresponding to adding a
handle). We seek to construct a family of solutions on $M_1\sharp M_2$,
such that for any compactly contained domains away from $p_1$ and
$p_2$, there are solutions in the family arbitrarily close to the
original data. The approximate solutions are obtained using conformal
blowups of $M_1\setminus \{p_1\}$ and $M_2\setminus \{p_2\}$, producing
solutions with asymptotically cylindrical ends. The connected sum is
then performed by cutting off far down on each cylinder and making
appropriate identifications. We assume that $K_j$ is nontrivial, and we
also make the nondegeneracy assumption that there are no nontrivial
conformal Killing vector fields on $M_j$ that vanish at~$p_j$. 

In normal coordinate neighborhoods $U_j\ni p_j$ around each point, the
metric takes the form $\gamma_j= \mathrm{d}r_j^2 + r_j^2 h_j(r_j)$,
where $h_j$ gives a family of metrics on $\mathbb{S}^2$, with
$h_j(0)=\r{g}_{\mathbb{S}^2}$ the unit round metric. For $j=1,2$, let
$\psi_j$ be a positive smooth function which is identically one on most
of~$M_j$, and which is $r_j^{1/2}$ near each $p_j$, say for
\xsmath{r_j< \r{r}}, and with $\psi_j=1$ outside of \xsmath{r_j<
2\r{r}} (take \xsmath{\mathrm{log}(\r{r}^{-2})>1}). With the change in
radial coordinate $t_j= -\log r_j$, we find $\psi_j^{-4} \gamma_j =
\mathrm{d}t_j^2 + \r{g}_{\mathbb{S}^2} +
O(\mathrm{e}^{-2t_j})$~\cite{cem}. One introduces a large parameter
$T\gg -2 \log$ \xsmath{\r{r}>1} to specify how far along each neck we
go before transitioning smoothly to an {\it exactly} cylindrical metric
$\r{\gamma} =\mathrm{d}s_j^2 + \r{g}_{\mathbb{S}^2}$. Indeed with
coordinate \xsmath{s_j= -\mathrm{log} (r_j/\r{r})- ({T}/{2})}, we have
$\psi_j^{-4} \gamma_j = \mathrm{d}s_j^2 + \r{g}_{\mathbb{S}^2}+
\mathrm{e}^{-T} \mathrm{e}^{-2s_j} \r{r}^2 \hat{h}_j$, and we can
smoothly transition in the region $(-1, -1/2) \times \mathbb{S}^2$ to
the exactly cylindrical metric on $(-1/2, \infty) \times \mathbb{S}^2$.
We appropriately form a quotient along the cylindrical pieces $(-1/2,
1/2) \times \mathbb{S}^2$ where the change of coordinates identifies
$s_1$ with $-s_2$; note that $s_j=0$ where \xsmath{r_j=\r{r}
\mathrm{e}^{-T/2}}, while $s_j= -{T}/{2}$ where \xsmath{r_j= \r{r}}.
With $s=s_1$ for $s\leq 0$ and $s=-s_2$ for $s\geq 0$, then on
$C_T\cong [-{T}/{2}, {T}/{2}] \times \mathbb{S}^2$, the metric is
$\gamma_T= \mathrm{d}s^2 + \r{g}_{\mathbb{S}^2}+ \mathrm{e}^{-T} \cosh
(2s) \hat{h}_T$, where $\hat{h}_T$ and its covariant derivatives with
respect to the cylindrical metric are bounded.

We have now put a metric on the connected sum, patching together the
conformal rescalings. Observe that $\psi^2_j \mu_j$ is TT for
$\psi_j^{-4} \gamma$, and just like for $\gamma_T$, we can patch
together the conformally rescaled TT-tensors along the cylinder, to
obtain an approximate solution~$\mu_T$. We want to conformally rescale
appropriately to get a suitable approximate solution to the constraint
equations. To do this, we patch the conformal factors $\psi_T :=
\chi_{1,T} \psi_1+ \chi_{2,T} \psi_2$, in such a way that $\psi_T$ is
the sum $\psi_1+\psi_2$ inside the region corresponding to $[-({T}/{2})
+ 1, ({T}/{2}) - 1] \times \mathbb{S}^2$, i.e.\ most of~$C_T$. In this
region, then
{$$
\psi_T= r_1^{1/2} + r_2^{1/2} =(\r{r}
\mathrm{e}^{-s-\frac{T}{2}})^{1/2} +(\r{r}
\mathrm{e}^{s-\frac{T}{2}})^{1/2} = 2 \r{r}^{1/2}
\mathrm{e}^{-\frac{T}{4}} \cosh\left(\frac{s}{2}\right).
$$}\unskip
In the exactly cylindrical piece $-1/2 < s < 1/2$, $\psi_T^4 \gamma_T$
is readily seen to have scalar curvature {\it zero}. By rotational
symmetry, then, it is {\it exactly Schwarzschild} in this region, with
mass \xsmath{2\r{r} \mathrm{e}^{-T/2}}. Thus for $T$ large, we have
a connected sum with the middle of the neck being a neighborhood of the
very small horizon sphere. 

From here, one can first modify $\mu_T$ to a $\gamma_T$-TT tensor
$\tilde{\mu}_T= \mu_T-\sigma_T$, with the perturbation $\sigma_T$ which
is globally small, with norm exponentially decaying in $T$
\cite[Proposition~5]{imp}. This step employs the conformal Killing
operator, and this is where the nondegeneracy assumption is used. The
proof is interesting and nontrivial, but we focus on the Lichnerowicz
operator here for the sake of brevity, noting its role in constant
scalar curvature gluing (see, e.g.~\cite{cem}).

Having sorted out $\tilde{\mu}_T$, one then solves the Lichnerowicz
equation
{$$
\mathcal{N}_T (u):=\Delta_{\gamma_T} u - \tfrac18 R(\gamma_T) u +
\tfrac18 |\tilde{\mu}_T|_{\gamma_T}^2 u^{-7} - \tfrac{1}{12}
\tau^2 u^5=0.
$$}\unskip
We make some remarks on the proof. 

First, $\psi_T$ can be used as an approximate solution, since one can
bound $\|\mathcal{N}_T(\psi_T)\|_{C^{k, \alpha}}\leq
C\mathrm{e}^{-T/2}$, where $C$ is independent of $T$ large
\cite[Proposition~6]{imp}, cf.~\cite[Proposition~4.6]{cem}. Indeed
outside of $C_T$, $\mathcal{N}_T(\psi_T)=0$. Now, where $\psi_T=
\psi_1+\psi_2$, we have \xsmath{\Delta_{\r{\gamma}} \psi_T - (1/8)
R(\r{\gamma}) \psi_T = 0} (we note \xsmath{R(\r{\gamma}) = 2}
for the cylindrical metric). As $\gamma_T$ is approximately cylindrical
in $[-1,1]\times \mathbb{S}^2$, with precise error estimates, and with
corresponding estimates on $\tilde{\mu}_T$ and $\psi_T$, we can
estimate $\mathcal{N}_T(\psi_T)$ in this region. The other regions in
$C_T$ are handled likewise, using that in $[-{T}/{2}, -1]\times
\mathbb{S}^2$ we have $\mathcal{N}_T(\psi_1)=0$, and in $[1, {T}/{2}]
\times \mathbb{S}^2$ we have $\mathcal{N}_T(\psi_2)=0$.

The next step is to get uniform (in $T$) bounds on the inverse of the
linearization $\mathcal{L}_T=D\mathcal{N}_T|_{\psi_T}$. We observe 
{\begin{equation} \label{eq:lich-lin}
\mathcal{L}_T = \Delta_{\gamma_T} - \tfrac18\left(R(\gamma_T)+ 7
|\tilde{\mu}_T|^2_{\gamma_T} \psi_T^{-8} + \tfrac{10}{3}
\tau^2\psi_T^4\right).
\end{equation}}\unskip
Showing that $\mathcal{L}_T$ is invertible for all $T$ large and
bounding the inverse are both handled by blowup arguments. Indeed,
suppose the existence of a sequence $T_m\nearrow \infty$, with $\sup_M
|\eta_m| =1$ but $\mathcal{L}_{T_m} \eta_m=0$. By the Schauder estimate
and reindexing a subsequence, we can assume that $\eta_m$ converges to
$\eta$ in~$C^2$. There are two basic scenarios: either we can also
arrange the subsequence so that $\sup |\eta_m|\geq c>0$ on
$M_1\setminus \{p_1\}$ or $M_2\setminus \{p_2\}$ (note that $(M_1\sharp
M_2, \gamma_T)$ contains more of $(M_1, \gamma_1) \cup (M_2, \gamma_2)$
isometrically as $T$ grows), or $\eta_m$ converges locally uniformly to
zero on these two pieces and the point where $\sup_M |\eta_m|=1$ is
attained travels further down the cylinder. By translation we can
center the supremum point, and obtain a nontrivial bounded $\eta$ on
the cylinder $\mathbb{R}\times \mathbb{S}^2$, satisfying
\xsmath{\Delta_{\r{\gamma}} \eta - ({1}/{4}) \eta =0}; this follows
from (\ref{eq:lich-lin}) by the estimates of $\tilde{\mu}_T$
and~$\psi_T$. As there are no nontrivial such $\eta$ (by separation of
variables or the maximum principle), we have a contradiction.

Thus we must be in the former case, say $\eta$ is nontrivial on
$M_1\setminus \{p_1\}$. In this case, $\gamma_T$ converges to
$\psi_1^{-4} \gamma_1$, $|\tilde{\mu}_T|^2_{\gamma_T}$ converges to
$\psi_1^{12}| \mu_1|^2_{\gamma_1}$, so that $\mathcal{L}_T$
(\ref{eq:lich-lin}) converges to
{$$
\tilde{\mathcal{L}}= \Delta_{\psi_1^{-4} \gamma_1} - \tfrac18
\left(R(\psi_1^{-4} \gamma_1) + 7 |\mu_1|^2_{\gamma_1} \psi_1^4 +
\tfrac{10}{3} \tau^2 \psi_1^4\right).
$$}\unskip
The conformal Laplacian satisfies a covariance property,
{$$
\Delta_{\psi_1^{-4} \gamma_1} \eta - \tfrac18 R(\psi_1^{-4} \gamma_1)
\eta = \psi_1^5 \left(\Delta_{\gamma_1} (\psi_1^{-1}\eta) - \tfrac18
R(\gamma_1) \psi_1^{-1} \eta\right),
$$}\unskip
while since $(\gamma_1, K_1= \mu_1+ ({\tau}/{3}) \gamma_1)$ solves the
vacuum constraints, we have $R(\gamma_1) = |\mu_1|_{\gamma_1}^2 -(2/3)
\tau^2$, and so we see that $\tilde{\mathcal{L}}\eta =0$ can be
re-written
{$$
\left(\Delta_{\gamma_1}- |\mu_1|_{\gamma_1}^2 - \tfrac13 \tau^2\right)
(\psi_1^{-1}\eta)=0.
$$}\unskip
The operator $(\Delta_{\gamma_1}- |\mu_1|_{\gamma_1}^2 - (1/3) \tau^2)$
has trivial kernel for nontrivial~$K_1$, which we assume. Thus it
admits a positive Green's function $G$ which has a pole or order
$r_1^{-1}$ at $p_1$, with $r_1$ the distance to~$p_1$. Since
$|\eta|\leq 1$ and $|\psi^{-1}| \leq r_1^{-1/2}$, $\psi^{-1} \eta$ must
extend smoothly across~$p_1$, and so we conclude it must vanish,
contradicting that $\eta$ is nontrivial on $M_1\setminus \{p_1\}$. 

Similar arguments can be used to control the inverse
\cite[Proposition~8]{imp} (cf.~\cite{stav-ppl, cem}). Writing a Taylor
expansion $\mathcal{N}_T(\psi_T+\eta) = \mathcal{N}_T (\psi_T) +
\mathcal{L}_T(\eta) + \mathcal{Q}_T(\eta)$, one can then show that for
large~$T$, the following map is a contraction on a suitable ball around
0 with small $T$-dependent radius in a H\"older space:
{$$
\eta \mapsto -\mathcal{L}_T^{-1} (\mathcal{N}_T(\psi_T) + \mathcal{Q}_T
(\eta)) = -\mathcal{L}_T^{-1} (N_T(\psi_T + \eta)- \mathcal{L}_T(\eta)).
$$}\unskip
The keys are that $\psi_T$ is a sufficiently good approximate solution,
and there is suitable control on $\mathcal{L}_T^{-1}$ uniform in $T$
large. 

The fixed point $\eta_T$ furnishes the desired conformal factor
$\tilde{\psi}_T=\psi_T+\eta_T$, and the solution $(\tilde{\psi}_T^4
\gamma_T, \tilde{\psi}_T^{-2} \tilde{\mu}_T + ({\tau}/{3})
\tilde{\psi}_T^4 \gamma_T)$ of the constraints on $M_1\sharp M_2$. The
difference between the solution and the starting data decays
exponentially in $T$ on compact subsets of each $M_j\setminus \{p_j\}$.

This result gives us a way to combine solutions to the vacuum
constraints (or to add a handle to a solution), with the constituent
parts approximately represented in the final solution. There are
analogous results in the asymptotically Euclidean and asymptotically
hyperboloidal setting, and an interesting
implementation~\cite{stav-ppl} of conformal gluing (including a
different{\break} family of approximate solutions) to construct initial
data with some features of point particle solutions. 

A natural question to ask is whether and under what conditions one can
preserve compact subdomains in the original data. Results in this
direction, see e.g.~\cite{cd, cip, cem}, use localized deformations, to
which we now turn. 

\section{Localized gluing} \label{sec:loc-glue}

In this section we will outline steps to glue an asymptotically flat
initial data set across a large annulus to an initial data set suitably
chosen from a given family of initial data. The perturbation will be
localized to the annular region. 

To bring the ideas across, we will present a proof sketch of the
following. For any $R>0$, let $\mathcal{E}_R=\{x\in \mathbb{R}^3:
|x|>R\}$, and let $\mathcal{E}= \mathcal{E}_1$. 

\begin{theorem} \label{th:em-ts}
Let $(g, E)$ be asymptoticaly flat, time-symmetric source-free initial
data for the Einstein--Maxwell equations on~$\mathcal{E}$, with ADM
mass $m(g)\neq 0$. For sufficiently large~$R$, there are source-free
initial data $(\bar{g}, \bar{E})$ on $\mathcal{E}$ which agree with the
given data data on $\mathcal{E} \setminus \mathcal{E}_R$, and on
$\mathcal{E}_{2R}$ agree with asymptotically flat data from a spacelike
slice in a suitable Reissner--Nordstr\"om spacetime, with the same
electric charge. 
\end{theorem}

The statement is for an asymptotic end (which could be outside
nontrivial charge distributions), whereas if $(M, g, E)$ comprised a
nontrivial complete asymptotically flat solution of the time-symmetric
Einstein--Maxwell constraints, then since $R(g) = 2\kappa |E|^2_g\geq
0$, we would have $m(g)>0$, by the Positive Mass Theorem. 

The corresponding theorem for the non-time symmetric case requires the
exterior slice to be chosen suitably in the Kerr--Newman family. For
simplicity of exposition, we focus on the time-symmetric case here. 
Before we begin, we note that non-vacuum gluing results using the
conformal method have been explored in~\cite{IMaxP}.

We remark that the argument will show the construction works for {\it
interior} gluing, i.e., the asymptotic exterior could be taken from the
given data, while the interior would then be from a suitable member of
the model family. For example, in the time-symmetric vacuum case, the
interior can be made to be a member of the Schwarzschild family. If
this interior solution has positive mass, then the resulting data will
be complete with two asymptotically flat ends and will contain the
Schwarzschild minimal sphere; as the gluing region will be far out in
the asymptotic end, this will be an outermost minimal surface in the
end, so the Penrose inequality restricts the mass based on that of the
exterior. In any case, we note that interior gluing has been used for
long-time existence results for the evolution problem~\cite{bieri-ptc,
ptc-LTIG}. 

\subsection{Proof of the Theorem} \label{sec3.1}

The method of proof we follow here is the same as in~\cite{cor:schw,
cs:ak}, cf.~\cite{cd}. We smoothly patch together the given data to
data from a slice in the model spacetime across an annular region
$\mathcal{A}_R=\{x: R< |x|<2R\}$ to form an approximate solution. As
$R$ increases, the metric approaches~\xsmath{\r{g}}, where
\xsmath{\r{g}_{ij} = \delta_{ij}} in the asymptotic coordinates, and
$E$ tends to zero (in the general case, $K$ also tends to zero). As
such, while the approximation improves with increasing~$R$, it
approaches a solution which admits nontrivial kernel of the formal
adjoint of linearization of the constraints. This finite-dimensional
kernel, arising from symmetries, is precisely the obstruction to a
local (and localized) perturbation result, as we now outline. 

The first key to achieving localized perturbation in our framework is a
weighted elliptic estimate (\ref{bewts}), which in the time-symmetric
case derives from an absolute estimate (\ref{bets2}) without boundary
terms. The estimates are used to guarantee that the solution obtained
lies in an appropriate weighted space, so that it extends smoothly
across~$\mathcal{A}_R$. Like the conformal method above, we need to get
uniform estimates on a family of linearized operators, but this can
only be done transverse to the finite-dimensional approximate kernel of
the adjoint of linearized operator (see estimate (\ref{cewts}) below). 
The presence of this approximate cokernel for large $R$ leads us to
employ a Lyapunov--Schmidt reduction: we use the weighted linear
elliptic estimates and a Picard-type quasi-Newtonian iteration to solve
a nonlinear projected problem that puts the constraint data into a
finite-dimensional space. We handle the obstruction space by choosing
a suitable member of the model family, in this case
Reissner--Nordstr\"om. 

To be precise, we recall the following operator introduced earlier:
{$$
\Psi(g, E) = (R(g) - 2 \kappa |E|_g^2, \mathrm{div}_g E).
$$}\unskip
For a perturbative analysis, we compute its linearization and formal
$L^2$-adjoint, for which one readily finds (and
recall~(\ref{eq:scal-lin}))
{$$
D\Psi_{(g, E)} (h, X) = \left(DR_g(h)- 4\kappa g(E, X) -2\kappa
h(E, E), \mathrm{div}_g(X)+\tfrac12 d(\mathrm{tr}_g h)(E)\right),
$$}\unskip
and hence 
{\begin{equation} \label{eq:lad-em}
D\Psi_{(g, E)}^*(f, \psi) = \left(DR_g^*(f) -2 \kappa f E^\flat \otimes
E^\flat- \tfrac12 \mathrm{div}_g (\psi E) g, -\mathrm{grad}_g \psi
- 4\kappa f E\right).
\end{equation}}\unskip

These simplify at the flat data, about which we will be perturbing. Let
\xsmath{\r{g}} be the Euclidean metric on $\mathbb{R}^3$, with
connection \xsmath{\r{\nabla}}. If we let \xsmath{D\r{\Psi} =
D\Psi_{(\r{g}, 0)}} with formal $L^2$-adjoint \xsmath{D \r{\Psi}^*},
then
{$$
D\r{\Psi} (h, X) = (DR_{\r{g}} (h),
\mathrm{div}_{\r{g}} X), \quad D\r{\Psi}^* (f, \psi)=
(DR^*_{\r{g}} (f), -\mathrm{grad}_{\r{g}} \psi).
$$}\unskip

The obstruction space is given by the elements $(f, \psi)$ in the
kernel of \xsmath{D\r{\Psi}^*}. As observed earlier, then $f$ must be a
constant plus a linear combination of the Cartesian coordinate
functions $x^1, x^2, x^3$, while $\psi$ must be a constant. The end of
the proof will show how to cover these directions by picking
appropriately the parameters (mass, center of mass and charge) of the
Reissner--Nordstr\"om data to which we glue, as we will detail in
Section~\ref{sec:cok}.

{\vspace*{2pt}}

\subsubsection{Weighted spaces} \label{sec3.1.1}

{\vspace*{2pt}}

We suppose $(\overline{\Omega}, \r{g})$ is compact, connected
with smooth nonempty boundary $\partial \Omega$. It suffices for our
purpose here to consider $\Omega$ to be an annular region. Given $N>0$,
and $0< r_1<r_0$, we let $\tilde{\rho}: \mathbb{R}\rightarrow
\mathbb{R}$ be a smooth nondecreasing function with $\tilde{\rho}(t)=0$
for $t\leq 0$, $\tilde{\rho}(t) = \mathrm{e}^{-N/t}$ for $0<t<r_1$, and
$\tilde{\rho}(t) = 1$ for $t\geq r_0$. We define $\rho$ on $\Omega$
by $\rho (x) = \tilde{\rho}(\mathrm{d}(x, \partial \Omega))$, where
$\mathrm{d}(x, \partial \Omega)$ is the \xsmath{\r{g}}-distance of $x$
to $\partial \Omega$; we choose $r_0$ suitably small, so that in
particular $\rho$ is smooth. On an annular region, $\rho$ decays to
zero at the boundary spheres, and is identically 1 near the middle of\break
the annulus.

Let $\|u\|^2_{L^2_{\rho}(\Omega)} = \int_{\Omega} |u|^2\,\rho\,
\mathrm{d}\mu_{\r{g}}$, while for $k$ a nonnegative integer, let
$\|u\|_{H^k_{\rho} (\Omega)}^2 = \sum_{\ell= 0}^k\|\r{\nabla}^{\ell}
u\|_{L^2_{\rho}(\Omega)}$. The resulting Hilbert spaces can in fact
can be identified as the closure of
$C^{\infty}(\overline{\Omega})$ in the relevant norm
\cite[Lemma~2.1]{cs:ak}. Furthermore by
\cite[Proposition~2.10]{cor-lan}, upon choosing $N$ suitably large,
there is a constant $C>0$ such that for all $u\in H^k_{\rho}(\Omega)$,
{\vspace*{2pt}}
{\begin{equation} \label{eqvwts}
\|u\rho^{1/2}\|_{H^k(\Omega)}\leq C \|u\|_{H^k_{\rho}(\Omega)}.
\end{equation}}\unskip
As such, for a bounded sequence in $H^k_{\rho}(\Omega)$, there is a
$u\in H^k_{\rho}(\Omega)$ and a subsequence $u_i$ for which $u_i
\rho^{1/2}$ converges to $u\rho^{1/2}$ in $H^{k-1}(\Omega)$ and weakly
in $H^k(\Omega)$, while $u_{i}$ converges weakly to $u$ in
$H^k_{\rho}(\Omega)$. We note that for metrics $g$ in a
$C^k(\overline{\Omega})$-bounded set, the norms defined as above but
using the measures $\mathrm{d}\mu_g$ and connection $\nabla_g$ are
uniformly equivalent.

The linearized operators enjoy estimates (\ref{bets2})--(\ref{bewts})
in these weighted spaces, which are used in solving the linearized
equation variationally (see Section~\ref{sec:pp}). There are also
weighted H\"older spaces built with the appropriate scaling to capture
the interior Schauder estimates (\ref{schwt})--(\ref{schwtinj}) that
follow for the solutions of the linearized equations. These spaces are
defined with a natural scaling compatible with the weight
function~$\rho$, and in fact there is a suite of weighted spaces
compatible with different weights~$\rho$. In~\cite{cor:schw}, we first
dealt with power decay, so $\rho = d^N$ near $\partial \Omega$, which
has a simple weighted estimate. While the exponential weighting was
discussed there, the weighted H\"older spaces used by Chru\'sciel and
Delay~\cite{cd} give a more elegant formulation of the interior
estimates, and allow for various kinds of weights all in one framework.

We follow an equivalent implementation as in~\cite{cor-lan}. We take a
smooth function $0<\phi<1$ with $\phi= d^2$ near the boundary, and with
$\overline{B_{\phi(x)}(x)} \subset \Omega$ for all $x\in \Omega$. We
let $\varphi$ be of the form $\varphi= \phi^r \rho^s$, with $r$ and $s$
real numbers, and let
{\vspace*{2pt}}
{$$
\|u\|_{C^{k, \alpha}_{\phi, \varphi}(\Omega)} = \sup_{x\in \Omega}
\left(\sum_{j=0}^k \varphi(x) \phi^j(x) \|\r{\nabla}^j
u\|_{C^0(B_{\phi(x)}(x))}+ \varphi(x) \phi^{k+\alpha}(x)
[\r{\nabla}^k u]_{0, \alpha; B_{\phi(x)}(x)}\right).
$$}\unskip
Every derivative is matched with a power of $\phi$, corresponding to
$\tilde{\rho}'(t) = N t^{-2} \tilde{\rho}(t)$, and by iterating this we
have $|\phi^k \rho^{-1} \nabla^k \rho|$ is bounded.

{\vspace*{2pt}}
\subsubsection{Integral estimates} \label{sec3.1.2}
{\vspace*{2pt}}

The key to time-symmetric estimates is that $DR_g^*$ is
overdetermined-elliptic. In fact
{\vspace*{2pt}}
{$$
\mathrm{Hess}_g f = DR_g^*(f) - \frac{1}{n-1} (\mathrm{tr}_g(DR_g^*
(f))+ fR(g)) g+ f \mathrm{Ric}(g),
$$}\unskip
from which we immediately see that there is a constant $C>0$ such that
on any open set $\Omega$, $\|f\|_{H^2(\Omega)} \leq C (\|D R_g^* (f,
\psi)\|_{L^2(\Omega)} + \|f\|_{H^1(\Omega)})$. 

With an application of Rellich's lemma we can replace the right-most
norm by $\|f\|_{L^2(\Omega)}$, and together with (\ref{eq:lad-em}) we
obtain
{\vspace*{2pt}}
{\begin{equation} \label{bets2}
\|f\|_{H^2(\Omega)} + \|\psi\|_{H^1(\Omega)} \leq C (\|D \Psi_{(g,
E)}^* (f, \psi)\|_{L^2(\Omega)}+\|f\|_{L^2(\Omega)}+\|\psi\|_{L^2
(\Omega)}).
\end{equation}}\unskip
 
We can replace the above norms with the weighted Sobolev norms by a
method used in~\cite{cor:schw}, namely we can replace $\Omega$ with
$\Omega_{\varepsilon}= \{x\in \Omega : d(x, \partial \Omega)>
\varepsilon\}$ for small $\varepsilon>0$, multiply the (squared)
estimates by $\tilde{\rho}'(\varepsilon)\geq 0$ and integrate: 
{\vspace*{2pt}}
{$$
\int_0^{\varepsilon_0} \tilde{\rho}'(\varepsilon) (\|f\|_{H^2
(\Omega_{\varepsilon})}^2 + \|\psi\|^2_{H^1(\Omega_{\varepsilon})})\,
\mathrm{d}\varepsilon \leq C \int_0^{\varepsilon_0} \tilde{\rho}'
(\varepsilon) (\|D\Psi_{(g, E)}^* (f, \psi)\|^2_{L^2(\Omega_{\varepsilon})}
+ \|f\|_{L^2(\Omega_{\varepsilon})}^2 + \|\psi\|^2_{L^2
(\Omega_{\varepsilon})})\,\mathrm{d}\varepsilon.
$$}\unskip
Integration by parts together with the coarea formula lead to 
{\vspace*{2pt}}
{\begin{equation}
\|(f, \psi)\|_{H^{2, 1}_{\rho}(\Omega)} := \|f\|_{H^2_{\rho}(\Omega)}+
\|\psi\|_{H^1_{\rho}(\Omega)} \leq C (\|D \Psi_{(g, E)}^*(f, \psi)
\|_{L^2_{\rho}(\Omega)} + \|f\|_{L^2_{\rho}(\Omega)} +
\|\psi\|_{L^2_{\rho}(\Omega)}). \label{bewts}
\end{equation}}\unskip
We note that (\ref{bets2}), and hence (\ref{bewts}), holds for {\it
any} data $(g, E)$, with a constant uniform for $(g, E)$ in a bounded
$C^2(\overline{\Omega})\times C^1(\overline{\Omega})$ neighborhood. 

While we focus on the time-symmetric case, we do note that the estimate
{\vspace*{2pt}}
{$$
\|Z\|_{H^1_{\rho}(\Omega)} \leq C (\|L_Z g\|_{L^2_{\rho}(\Omega)} +
\|Z\|_{L^2_{\rho}(\Omega)})
$$}\unskip
for $Z$ is more subtle. The proof from \cite[Lemma~5.1]{cor-lan} is a
variation on the argument from \cite[Lemma~4.1]{cs:ak}; both use the
elementary inequality \xsmath{\Delta_{\r{g}} \rho \geq (1/2) N^2 d^{-4}
\rho} near $\partial \Omega$ and an integration-by-parts argument, and
the estimate holds for $N$ sufficiently large. The analogue of
(\ref{bewts}) for the full constraints operator follows.

\subsubsection{Uniform estimate transverse to cokernel}
\label{sec3.1.3}

We can promote (\ref{bewts}) to a coercivity estimate, at least in
directions transverse to the cokernel $\mathcal{K}$ at flat data.
Recall $\mathcal{K}= \mathrm{span}\{1, x^1, x^2, x^3\} \oplus
\mathrm{span}\{1\} \subset H^{2}_{\rho}(\Omega)\times
H^{1}_{\rho}(\Omega) =: H^{2, 1}_{\rho}(\Omega)$. Fix a nontrivial
smooth bump function $\zeta$ supported on~$\Omega$. Let $S_g$ be the
$L^2(\Omega, \mathrm{d}\mu_g)$-orthogonal complement of~$\zeta
\mathcal{K}$.
 
We claim there is a $C>0$ and a $C^2(\overline{\Omega}) \times
C^1 (\overline{\Omega})$ neighborhood $\mathcal{U}$ of the flat data
\xsmath{(\r{g}, 0)} such that for all data $(g, E)\in \mathcal{U}$ and
all $(f, \psi)\in S_g$,
{\begin{equation} \label{cewts}
\|(f, \psi)\|_{H^{2, 1}_{\rho}(\Omega)} \leq C \|D\Psi_{(g, E)}^*(f,
\psi)\|_{ L^2_{\rho}(\Omega)}.
\end{equation}}\unskip
The same holds if $S_g$ is replaced by any closed subspace $S$ of
$H^{2, 1}_{\rho}(\Omega)$ with $S\cap \mathcal{K}=\{0\}$. 

There are a couple ways one can prove this. For instance, one can argue
along the lines of the proof of \cite[Proposition~3.1]{cem},
establishing the {\it unweighted} analogue of (\ref{cewts}) on
$\Omega_{\varepsilon}$, with a uniform constant for $\varepsilon$
small, and proceed as in the proof of (\ref{bewts}). Another approach
is to use (\ref{bewts}) along with a standard compactness argument:
from a sequence $(f_i, \psi_i)\in S_{g_i}$, with \xsmath{(g_i, E_i)\rightarrow
(\r{g}, 0)}, $\|(f_i, \psi_i)\|_{H^{2, 1}_{\rho}(\Omega)} =1$
and $\|D\Psi_{(g, E)}^*(f_i, \psi_i)\|_{L^2_{\rho}(\Omega)}\rightarrow
0$, one can construct a nontrivial $(f, \psi)\in S_{\r{g}} \cap
\mathcal{K}$, which furnishes a contradiction. Indeed as noted after
(\ref{eqvwts}), there is a $(f, \psi)\in H^{2, 1}_{\rho}(\Omega)$ such
that, upon re-indexing a subsequence, $(f_i, \psi_i)\rho^{1/2}$
converges to $(f, \psi) \rho^{1/2}$ in $H^{1}(\Omega)\times
L^2(\Omega)$. In particular, $(f_i, \psi_i)$ converges in
$L^2_{\mathrm{loc}}(\Omega)$ to $(f, \psi)$, so that
\xsmath{D\r{\Psi}^* (f, \psi)= 0} weakly. Hence $(f, \psi)\in
\mathcal{K}$. Since $(f_i, \psi_i)$ is
$L^2(\mathrm{d}\mu_{g_i})$-orthogonal to $\zeta \mathcal{K}, (f,
\psi) \in S_{\r{g}}$. Moreover, since $(f_i, \psi_i)\rightarrow
(f, \psi)$ in $ L^2_{\rho}(\Omega)$, we see from (\ref{bewts}) that
$(f, \psi)$ is nontrivial, which is a contradiction.

\subsubsection{The projected problem} \label{sec:pp}

With $\zeta$, $\mathcal{K}$ and $S_g$ as above, let $\Pi_g$ be the
$L^2(\Omega, \mathrm{d}\mu_g)$-orthogonal projection 
to $S_g$, and let \xsmath{\r{\Pi} =
\Pi_{\r{g}}}. For simplicity, we may arrange that $\zeta$ is
supported where $\rho=1$. We study the map \xsmath{(g, E) \mapsto
\r{\Pi} (\Psi(g, E))}, with linearization \xsmath{(h, X) \mapsto
\r{\Pi} (D\Psi_{(g, E)}(h, X))}.

The goal is to solve $\Psi(g, E)=0$. We start by solving a projected
problem:
{$$
\r{\Pi} (\Psi(g_0+h, E_0+X))= 0,
$$}\unskip
where $(g_0, E_0)$ is close to \xsmath{(\r{g}, 0)}. The linearization
of this is 
{\begin{equation} \label{eq:prlin0}
\r{\Pi} (D\Psi_{(g_0, E_0)}(h, X))= -\r{\Pi} (\Psi(g_0,
E_0)).
\end{equation}}\unskip
To solve this we consider the functional on $H^{2,1}_{\rho}(\Omega)
\cap S_{g_0}$ given by
{$$
\mathcal{G}_0 (f, \psi) = \int_{\Omega} \left(\frac12 \rho |
D\Psi^*_{(g_0, E_0)} (f, \psi)|_{g_0}^2 + \Pi_{g_0} (\Psi(g_0,
E_0)) \cdot (f, \psi) \right) \mathrm{d}\mu_{g_0}.
$$}\unskip
If $\mathcal{G}_0$ is stationary at $(f_0, \psi_0)$, the
Euler--Lagrange equation
{$$
\Pi_{g_0} (D\Psi_{(g_0, E_0)} \rho D\Psi^*_{(g_0, E_0)}(f_0,
\psi_0))= -\Pi_{g_0} (\Psi(g_0, E_0))
$$}\unskip
is implied by the following, which holds for all $(u, v)$ compactly
supported:
{$$
\int_{\Omega} (\Pi_{g_0} (D\Psi_{(g_0, E_0)} \rho D\Psi^*_{(g_0, E_0)}
(f, \psi)) + \Pi_{g_0} (\Psi(g_0, E_0))) \cdot (u, v)\,\mathrm{d}
\mu_{g_0}.
$$}\unskip
This follows trivially for $(u,v)\in \zeta \mathcal{K}$, while for
$(u,v)\in S_{g_0}$, the leading projection operator in the above
integral equation can be removed, and so it follows by stationarity of
the functional at $(f, \psi)\in H^{2, 1}_{\rho}(\Omega) \cap S_{g_0}$.

We conclude that $\Psi(g_0, E_0)+ D\Psi_{(g_0, E_0)} \rho
D\Psi^*_{(g_0, E_0)}(f_0, \psi_0) \in \zeta \mathcal{K}$, and so we see
that (\ref{eq:prlin0}) is solved with $(h, X)= \rho D\Psi^*_{(g_0,
E_0)}(f_0, \psi_0)$. Moreover, combining (\ref{cewts}) with
$\mathcal{G}_0(f_0, \psi_0)\leq 0$ we find
{\begin{equation} \label{eq:varest}
\|(f_0, \psi_0)\|_{H^{2, 1}_{\rho}(\Omega)} \leq C \|\Pi_{g_0}
(\Psi (g_0, E_0))\|_{L^2_{\rho^{-1}}(\Omega)}\leq C'\| \r{\Pi}
(\Psi (g_0, E_0))\|_{L^2_{\rho^{-1}}(\Omega)}.
\end{equation} }\unskip
 
\subsubsection{Pointwise estimates} \label{sec3.1.5}

The above shows we can solve the desired equation at the linear level,
with an integral estimate. To iterate linear corrections to solve the
nonlinear problem, we require a pointwise estimate of the variational
solution, in part to establish the convergence of a quasi-Newtonian
iteration scheme. 

The operator $(f, \psi) \mapsto D\Psi_{(g,E)} D\Psi^*_{(g, E)}(f,
\psi)$, and in fact for our purposes, the weighted version
{\begin{equation} \label{wtdop}
(f, \psi) \mapsto \rho^{-1} D\Psi_{(g,E)}\rho D\Psi^*_{(g, E)}
(f, \psi)=: (L_{11} f + L_{12}\psi, L_{21}f + L_{22} \psi),
\end{equation}}\unskip
is an elliptic operator of mixed order that admits Douglis--Nirenberg
weights~\cite{DN-INT}: $L_{ij}$ is order $s_i+t_j$, where $s_1= 0$,
$s_2=-1$, $t_1=4$, $t_2=3$, i.e.\ $L_{11}$ is order four, $L_{12}$ and
$L_{21}$ are order three, and $L_{22}$ is order two. For example, at
the flat data \xsmath{(\r{g}, 0)}, \xsmath{D\r{\Psi}
D\r{\Psi}^*(f, \psi) = (2\Delta_{\r{g}}^2 f, -
\Delta_{\r{g}} \psi)}. The diagonal terms of (\ref{wtdop}) have
the analogous principal parts, but there are also off-diagonal terms
(though they are small for $E$ small), and terms from the weight, which
can in principle blow up on approach to the boundary at a rate which is
a power of $d^{-1}$. This behavior can be accommodated in the interior
elliptic estimates by a scaling argument, which is efficiently coded by
the weighted H\"older norms introduced earlier. Indeed if we write
$L_{ij}$ in local coordinates as $L_{ij} = \sum_{|\beta|\leq s_i+t_j}
b_{ij}^\beta \partial^\beta$, then it is not hard to see that for $0<
\alpha<1$, there is a $C>0$ such that $\|b_{ij}^\beta\|_{C^{-s_i,
\alpha}_{\phi, \phi^{s_i+t_j-|\beta|}}(\Omega)}\leq C$. We record the
following interior Schauder estimate (see \cite[Theorem~5.8]{cor-lan},
e.g.), with $\varphi_j=\phi^{(3/2)+4-t_j}\rho^{1/2}$, i.e.\ $\varphi_1=
\phi^{3/2}\rho^{1/2}$ and $\varphi_2= \phi^{5/2}\rho^{1/2}$:
{\begin{eqnarray} 
&& \|f\|_{C^{4, \alpha}_{\phi, \varphi_1}(\Omega)} + \|\psi\|_{C^{3,
\alpha}_{\phi, \varphi_2}(\Omega)}\nonumber\\
&& \quad \leq C (\|L_{11} f + L_{12}\psi
\|_{C^{0, \alpha}_{\phi, \phi^4 \varphi_1} (\Omega)}+ \|L_{21} f+
L_{22}\psi \|_{C^{1, \alpha}_{\phi, \phi^2\varphi_2} (\Omega)}+
\|f\|_{L^2_{\rho}(\Omega)} + \|\psi\|_{L^2_{\phi^2\rho}(\Omega)}).
\label{schwt}
\end{eqnarray}}\unskip

We apply this to get a pointwise estimate of $(f_0, \psi_0)$, where
$(h, X)= \rho D\Psi^*_{(g_0, E_0)}(f_0, \psi_0)$ is the variational
solution to~(\ref{eq:prlin0}). The integral terms can be estimated
by~(\ref{eq:varest}). We can also replace terms in (\ref{schwt})
involving the operator by the their projected values: since $\rho=1$ on
the support of~$\zeta$, 
{$$
\rho^{-1} D\Psi_{(g_0, E_0)} \rho D\Psi^*_{(g_0, E_0)}(f_0, \psi_0)-
\rho^{-1} \r{\Pi} (D\Psi_{(g_0, E_0)} \rho D\Psi^*_{(g_0, E_0)}
(f_0, \psi_0))\in \zeta \mathcal{K}.
$$}\unskip
Using the fact that $\zeta \mathcal{K}$ is finite-dimensional (so all
norms on it are equivalent) and an elementary compactness argument, one
can derive (cf.\ the proof of \cite[Proposition~6.4]{cor-lan})
{\begin{eqnarray}
&& \|f_0\|_{ C^{4, \alpha}_{\phi, \varphi_1}(\Omega)} +
\|\psi_0\|_{C^{3, \alpha}_{\phi, \varphi_2}(\Omega)} +\|(f_0,
\psi_0)\|_{H^{2, 1}_{\rho}(\Omega)}\nonumber\\
&& \quad \leq C (\|\rho^{-1} \r{\Pi}(\Psi(g_0, E_0))
\|_{C^{0, \alpha}_{\phi, \phi^4\varphi_1} (\Omega)\times C^{1,
\alpha}_{\phi, \phi^2\varphi_2} (\Omega)} + \|f_0\|_{L^2_{\rho}
(\Omega)} + \|\psi_0\|_{L^2_{\phi^2\rho}(\Omega)})\nonumber\\
&& \quad \leq C (\| \r{\Pi}(\Psi(g_0, E_0)) \|_{C^{0,
\alpha}_{\phi, \phi^4\varphi_1 \rho^{-1}} (\Omega)\times C^{1,
\alpha}_{\phi, \phi^2\varphi_2 \rho^{-1}} (\Omega)} +
\|\r{\Pi} (\Psi (g_0, E_0))\|_{L^2_{\rho^{-1}}(\Omega)})
\label{schwtinj}
\end{eqnarray}}\unskip
where we used the variational estimate (\ref{eq:varest}). We define
norms such that we can abbreviate the above result as
{\begin{equation} \label{pvest}
\|(f_0, \psi_0)\|_{\mathcal{B}^{4, 3}} \leq C \|\r{\Pi}
(\Psi(g_0, E_0))\|_{\mathcal{B}^{0,1}}.
\end{equation}}\unskip

\subsubsection{Iteration} \label{sec3.1.6}

If we let $(h_0, X_0) = \rho D\Psi^*_{(g_0, E_0)}(f_0, \psi_0)$ and
define
{$$
\|(h, X) \|_{\mathcal{B}^{2,2}} = \|(h, X)\|_{C^{2, \alpha}_{\phi,
\phi^{\frac{7}{2}} \rho^{-\frac12}} (\Omega)} + \|(h, X)
\|_{L^2_{\rho^{-1}}(\Omega)},
$$}\unskip
it follows from (\ref{pvest}) that 
{$$
\|(h_0, X_0)\|_{\mathcal{B}^{2,2}} \leq C \|(f_0, \psi_0)
\|_{\mathcal{B}^{4, 3}}\leq C'\| \r{\Pi}(\Psi(g_0, E_0))
\|_{\mathcal{B}^{0,1}}.
$$}\unskip

The point is that by design, \xsmath{\r{\Pi}(\Psi(g_0+h_0, E_0+X_0))}
is quadratic in $(h_0,X_0)$, i.e.\ in \xsmath{\r{\Pi}(\Psi(g_0, E_0))}.
We set up a recursion by a quasi-Newtonian scheme to determine a
sequence $(h_m, X_m)=\rho D\Psi^*_{(g_0, E_0)}(f_m, \psi_m)$,
satisfying 
{$$
\r{\Pi} (D\Psi_{(g_0, E_0)} \rho D\Psi^*_{(g_0, E_0)}(f_{m+1},
\psi_{m+1})) = -\r{\Pi} (\Psi(g_m, E_m)),
$$}\unskip
where $g_m=g_0 + \sum_{k=0}^{m-1} h_k$, and $E_m= E_0 +
\sum_{k=0}^{m-1} X_k$. Note that the linearization is computed at the
starting data, which we take to be smooth. 

For sufficiently small $\Psi(g_0, E_0)$, the resulting convergence is
not quadratic, but it is super-linear; for more details, see
\cite[Section~6.3 and Appendix~D]{cor-lan}, cf.~\cite{cor:schw, cs:ak},
and \cite[Appendix~G]{cd}. We obtain $g_m\rightarrow g= g_0+h$ and
$E_m\rightarrow E = E_0 + X$, with \xsmath{\r{\Pi}(\Psi (g, E))= 0}, in
other words, $\Psi(g, E)\in \zeta \mathcal{K}$. Here $(h, X) = \rho
D\Psi^*_{(g_0, E_0)}(f, \psi)\in \mathcal{B}^{2,2}$, where $f=
\sum_{k=0}^\infty f_k$ and $\psi = \sum_{k=0}^\infty \psi_k$, with
$(f, \psi) \in \mathcal{B}^{4,3} \cap S_{g_0}$, and $\|(h,
X)\|_{\mathcal{B}^{2, 2}} \leq C \|\r{\Pi} (\Psi(g_0,
E_0))\|_{\mathcal{B}^{0,1}}$. By elliptic bootstrapping on the
quaslinear equation $\Psi(g, E)\in \zeta \mathcal{K}$, one can show
that $(h, X)$ extends smoothly by zero across $\partial \Omega$.

\subsubsection{Handling the cokernel} \label{sec:cok}

Thus far we have not glued anything, we have just referenced data near
\xsmath{(\r{g}, 0)} on a domain~$\Omega$. We will now restrict to
$\Omega= \mathcal{A}_1$, the unit annulus. Given asymptotically flat
initial data $(g, E)$ satisfying $\Psi(g, E) = 0$, we use a cutoff
function to smoothly patch $(g, E)$ to $(g_{\mathrm{RN}},
E_{\mathrm{RN}})$ in $\mathcal{A}_R$, obtaining $(\hat{g}_R,
\hat{E}_R)$. We may assume for exposition that $g_{ij}(x) -\delta_{ij}
= O (|x|^{-1})$ and $E^i(x)= O(|x|^{-2})$ (and derivatives used in the
proof attain additional decay), though we could formulate somewhat more
generally. For instance, we could pose the starting solution of the
constraints is asymptotically flat of order~$q$, i.e.\ in asymptotic
coordinates $|x|>1$, $\partial_x^\alpha (g_{ij}(x)-\delta_{ij}) = O
(|x|^{-q-|\alpha|})$ and $\partial_x^\beta E^i=
O_1(|x|^{-1-q-|\beta|})$ for $q>1/2$, $|\alpha|\leq 2$ and $|\beta|\leq
1$. Notice that $R(g) = 2\kappa |E|_g^2= O (|x|^{-2-2q})$ which is
integrable, and hence the ADM mass exists. In any case, we consider
$q=1$ for exposition, and the interested reader can readily generalize.

The exterior we attach depends on the parameters $m$, $Q$, as well as
$c\in \mathbb{R}^3$ about which we center $(g_{\mathrm{RN}}, E_{\mathrm{RN}})$. 
We then pullback to the unit annulus, as follows:
for $x \in\mathcal{A}_1$, let $(\tilde{g}_{R})_{ij}(x)=
(\hat{g}_R)_{ij} (Rx) $, and $\tilde{E}^i_R(x) = R\hat{E}^i_R(Rx)$. 
Observe that $(\tilde{g}_R, \tilde{E}_R)$ differs from \xsmath{(\r{g},
0)} in $\mathcal{A}_1$ by $O(R^{-1})$, while $\Psi(\tilde{g}_R,
\tilde{E}_R)(x) = R^2 \Psi(\hat{g}_R, \hat{E}_R)(Rx)= O(R^{-1})$. We
see the scaling was chosen so that we can solve the vacuum constraints
on $\mathcal{A}_1$ and simply scale back. Moreover, for all large~$R$,
we can apply the previous analysis to find $(h_R, X_R)$ that extend
smoothly by zero outside $\mathcal{A}_1$, such that with
$(\overline{g}_R, \overline{E}_R)= (\tilde{g}_R+ h_R, \tilde{E}_R+
X_R)$, we have $\r{\Pi}(\Psi(\overline{g}_R,\overline{E}_R))=0$.

The final step is to arrange $\Psi(\overline{g}_R, \overline{E}_R)=0$,
for which we recall $\mathcal{K} = \mathrm{span}\{1, x^1, x^2, x^3\}
\oplus \mathrm{span} \{1\}$. Projecting $\Psi(\overline{g}_R,
\overline{E}_R)$ onto the second summand yields (by scaling) 
{$$
R\int_{\mathcal{A}_1} \mathrm{div}_{\overline{g}_R} \overline{E}_R\,
\mathrm{d}\mu_{\overline{g}_R}=R \int_{\partial \mathcal{A}_1}
\tilde{E}_R \cdot \nu_{\tilde{g}_R}\,\mathrm{d}\sigma_{\tilde{g}_R}=
\int_{\{|x|=2R\}} E_{\mathrm{RN}}\cdot \nu_{g_{\mathrm{RN}}}\,
\mathrm{d}\sigma_{g_{\mathrm{RN}}} - \int_{\{|x|=R\}} E \cdot \nu_g\,
\mathrm{d}\sigma_g.
$$}\unskip
Since $\mathrm{div}_g E=0$, the integral $\int_{\{|x|=R\}} E \cdot
\nu_g\,\mathrm{d}\sigma_g$ is independent of $R>1$, and the value is
$4\uppi Q$. Thus if we patch the solution to Reissner--Nordstr\"om data
with charge~$Q$, then we conclude $\int_{\mathcal{A}_1}
\mathrm{div}_{\overline{g}_R} \overline{E}_R\,
\mathrm{d}\mu_{\overline{g}_R}=0$. Since
$\mathrm{div}_{\overline{g}_R} \overline{E}_R = b\zeta$ for a
constant~$b$, we conclude $b=0$, and hence
$\mathrm{div}_{\overline{g}_R} \overline{E}_R=0$. The reader might
compare this with our earlier remark in Section~\ref{sec:ece} about how
to prescribe the electric field data. 

For the rest of the section, we extend the summation convention to all
pairs of repeated indices (even if both are down). Since in Cartesian
coordinates, \xsmath{DR_{\r{g}} (h)= h_{ij, ij} - h_{ii, jj}}, we see
that by simple expansion of the scalar curvature term 
{\begin{eqnarray*} 
\int_{\mathcal{A}_1}(R(\overline{g}_R) - 2 \kappa
|\overline{E}_R|^2_{\overline{g}_R})\,\mathrm{d}x &=& \int_{\partial
\mathcal{A}_1} ((\tilde{g}_R)_{ij, i} - (\tilde{g}_R)_{ii, j})\,
\r{\nu}^j\,\mathrm{d}\r{\sigma} + O(R^{-2})\\
&=& R^{-1} \int_{\partial \mathcal{A}_R} ((\hat{g}_R)_{ij, i} -
(\hat{g}_R)_{ii, j})\,\r{\nu}^j\,\mathrm{d}\r{\sigma}+
O(R^{-2})\\
&=& R^{-1}16 \uppi (m - m(g)) + O(R^{-2}).
\end{eqnarray*}}\unskip
So the projection is to leading order controlled by the mass $m$ of the
Reissner--Nordstr\"om exterior. 

The projection onto the linear elements requires a bit more attention. 
For the metric~$g$, we find that for any $r>r_0>1$, 
{$$
\int_{\{r_0\leq |x| \leq r\}} x^k (g_{ij, ij} - g_{ii, jj})\,\mathrm{d}
x = \left(\int_{\{|x|=r\}} - \int_{\{|x|= r_0\}}\right)(x^k (g_{ij, i}-
g_{ii, j})- (\delta^{k}_{i} g_{ij} - \delta^k_{j} g_{ii}))
\r{\nu}^j\,\mathrm{d}\r{\sigma}.
$$}\unskip
Note that in the last terms of the integral on the right, we could
replace the terms of the form $g_{\ell m}$ with $g_{\ell m}-
\delta_{\ell m}$. Moreover, {\it using the constraint equations}, the
integrand $(g_{ij, ij} - g_{ii, jj})= O(|x|^{-4})$, and so the left
hand side is $O(\log r)$.

If instead the above is evaluated at the Reissner--Nordstr\"om data, by
using the asymptotic parity and the fact that the center of mass does
not appear in the $|x|^{-1}$ term in the expansion, we can show that
for $|c|\leq {r}/{2}$ and for $m$ bounded, $\int_{\{|x| \geq r\}} x^k
((g_{\mathrm{RN}})_{ij, ij} -(g_{\mathrm{RN}})_{ii, jj})\,\mathrm{d}x
= O(r^{-1})$, on the one hand, while on the other
{$$
\lim_{r\rightarrow \infty} \int_{\{|x|= r\}}(x^k ((g_{\mathrm{RN}})_{ij,
i}-(g_{\mathrm{RN}})_{ii, j})-(\delta^{k}_{i} (g_{\mathrm{RN}})_{ij} -
\delta^k_{j} (g_{\mathrm{RN}})_{ii})) \r{\nu}^j\,
\mathrm{d}\r{\sigma}=16\uppi m c^k.
$$}\unskip

Assembling all of this, noting that the rescaled exterior is a
Reissner--Nordstr\"om with center $\tilde{c} = cR^{-1}$, we obtain
{\begin{eqnarray}
&& \int_{\mathcal{A}_1} x^k(R(\overline{g}_R) - 2 \kappa
|\overline{E}_R|^2_{\overline{g}_R})\,\mathrm{d}x\nonumber\\
&& \quad = \int_{\partial \mathcal{A}_1} x^k ((\tilde{g}_R)_{ij, i}
- (\tilde{g}_R)_{ii, j})\, \r{\nu}^j\,\mathrm{d}\r{\sigma}
+ O(R^{-2})\nonumber\\
&& \quad = R^{-2} \int_{\{|x|=2R\}}(x^k((g_{\mathrm{RN}})_{ij, i}-
(g_{\mathrm{RN}})_{ii,j})- (\delta^{k}_{i} (g_{\mathrm{RN}})_{ij} -
\delta^k_{j}(g_{\mathrm{RN}})_{ii})) \r{\nu}^j\,\mathrm{d}
\r{\sigma}\nonumber\\
&& \qquad -\;R^{-2}\int_{\{|x|= R\}} (x^k(g_{ij, i}-g_{ii, j})-
(\delta^{k}_{i} g_{ij} - \delta^k_{j} g_{ii})) \r{\nu}^j\,
\mathrm{d}\r{\sigma} + O(R^{-2})\nonumber\\
&& \quad = R^{-1} \left(16\uppi m \tilde{c}^k + O\left(\frac{\log R}
{R}\right)\right) + O(R^{-2}). \label{eq:mc}
\end{eqnarray}}\unskip
Thus we get a map
{$$
(m, \tilde{c}) \mapsto \frac{R}{16\uppi} \int_{\mathcal{A}_1} (1, x)(R
(\overline{g}_R) - 2 \kappa |\overline{E}_R|^2_{\overline{g}_R})\,
\mathrm{d}x = (m-m(g) + \xi(m, \tilde{c}), m \tilde{c} +\Xi
(m, \tilde{c})),
$$}\unskip
where $\xi(m, \tilde{c}) = O(R^{-1})$ and $\Xi(m, \tilde{c})= O({\log
R}/{R})$. Moreover, the construction is {\it continuous} with respect
to $m$ and $\tilde{c}$, and so $\xi$ and $\Xi$ are continuous; this is
fairly straightforward but cumbersome to prove (see, e.g.\ 
\cite[Section~A.7]{acp:multi}). A fixed-point or degree argument can
finish the proof from here. For example, the function
{$$
F(m, \tilde{c}) = \left(m(g) - \xi(m, \tilde{c}), -\frac{\Xi(m,
\tilde{c})}{m(g) - \xi(m, \tilde{c})}\right),
$$}\unskip
is readily seen to map $\{(m, \tilde{c}) : |m-m(g)|\leq {|m(g)|}/{2},
|\tilde{c}|\leq 1/2\}$ to itself, and the Brouwer fixed point theorem
can then be applied. The fixed point gives parameters that solve the
problem. 

We remark that we did not assume the original metric $g$ enjoyed enough
parity symmetry to admit a well-defined center of mass. If this were
the case, we could argue that the resulting center $c= R\tilde{c}$ is
close to that of the original, cf.~\cite{cor:schw, cs:ak}. 

\subsection{Further applications of deformation} \label{sec3.2}

Suppose the starting data, say $(g_0, E_0)$, is such that $D\Psi_{(g_0,
E_0)}^*$ admits only trivial kernel in~$\Omega$. The same perturbative
analysis above can be applied, the only difference is that since there
is no kernel to deal with, we do not need the projected operator, and
so the technique yields a localized deformation result: given a
compactly contained subdomain of $\Omega$ and sufficiently small
compactly supported $S$ and~$\sigma$, say, there are $h$ and $X$ which
smoothly extend by zero across $\partial \Omega$ with $\Psi(g_0+h,
E_0+X)= \Psi(g_0, E_0) + (S, \sigma)$. Similar comments apply to
scalar curvature deformation~\cite{cor:schw}, as well as for the
general constraints case~\cite{cs:ak}, though in the latter case
maintaining the dominant energy condition at the borderline is somewhat
subtle~\cite{cor-lan}. 

As mentioned earlier, one can combine conformal gluing and localized
gluing to obtain connected sum constructions leaving the original data
intact outside the gluing neck~\cite{cd, cip, cem}. The additional
assumption needed, beyond any needed for the conformal gluing, is that
the corresponding operator for the localized gluing problem has trivial
kernel. If you have a sequence of solutions of the nonlinear problem
(vacuum constraints, or a scalar curvature equation, say) which
converges to a solution for which the linearized adjoint has trivial
kernel in~$\Omega$, then you can {\it glue back in} the original data
inside some appropriate subdomain, and use the localized perturbations
(unobstructed in this case) to reimpose the nonlinear equation. 

For an interesting application combining localized deformations of
constraint data with conformal deformations, consider an asymptotically
flat metric $g$ with nonnegative scalar curvature $R(g)$, which is not
{\it static} outside of some bounded region $\Omega\subset
\mathbb{R}^3$ (in the sense that the adjoint $DR_g^*$ of the
linearization of the scalar curvature has trivial kernel in the
exterior $\mathbb{R}^3\setminus \overline{\Omega}$). It is possible
using localized deformations to find a compactly contained domain
$V\subset R^3\setminus \overline{\Omega}$ such that for $S$ compactly
supported in~$V$, and for $\varepsilon>0$ sufficiently small, there is
a metric $\bar{g}_{\varepsilon}$ with $R(\bar{g}_{\varepsilon}) = R(g)
+ \varepsilon S$, and $\bar{g}_{\varepsilon}- g=O(\epsilon)$ is
supported in~$V$. If $S\geq 0$ and is nontrivial, one can then use an
elementary conformal transformation to {\it bring the ADM mass down}
while keeping the scalar curvature nonnegative, with conformal factor
$u_{\varepsilon}>0$ which is one in a neighborhood of
$\overline{\Omega}$~\cite{jc-nbm}. In terms of asymptotically flat
extensions of a given region, then, this shows that unless the data
outside the region were static in the above sense, there would be
extensions with nonnegative scalar curvature and lower ADM mass, and
hence in particular the original configuration does not minimize the
Bartnik mass. To summarize, a localized deformation is used to bump the
scalar curvature up (adding energy density) in a compact set outside
the region $\overline{\Omega}$, which does not change the ADM mass, and
then a conformal deformation removes some of this scalar curvature and
in so doing brings down the ADM mass, leaving the original geometry{\break} of
$\overline{\Omega}$ intact.

Similar considerations hold for the general (non-time-symmetric) case,
though maintaining the dominant energy condition seems to require a
modification of the constraints operator to have trivial
kernel~\cite{cor-lan, jc-snbm}. An interpretation of the kernel
elements for the modified operator, akin to that in~\cite{monc:stat},
and has been established in the recent work~\cite{huang-lee-bm} on
constraint deformations and the dominant energy condition, with
striking applications to the Bartnik mass. 

Another application is the gluing of initial data to {\it interpolate}
their scalar curvatures~\cite{delay-int}, or more generally
interpolating $(\mu, J)$ from the constraints operator~\cite{cor-lan}. 
As an interesting example, one can glue initial data for a Kerr
solution (including Schwarzschild) to a different Kerr solution across
an annular region where the dominant energy condition can be maintained
(see Remark after Theorem~1.4 in~\cite{cor-lan}).

Finally, we discussed earlier the use of the conformal method to
achieve solutions with certain prescribed asymptotic behavior, with
stability for the Einstein evolution in mind. Motivated by questions
about the asymptotics posed in~\cite{car-sch}, the work~\cite{lef-ngu}
proceeds to seed the desired asymptotics, and then correct to solutions
of the constraints, but in a different manner than that
of~\cite{FST-mink, FST-BH}, in part employing variational methods of
the kind used in~\cite{cs:ak, car-sch}. 

\subsection{Alternate approach to solving the projected problem}
\label{sec3.3}

We review the steps in the procedure shown above for gluing to an
asymptotically flat model family. We start by patching initial data to
a member of the model family, with the following parameters at our
disposal: $R$ (measuring the scale), the exterior mass $m$ and center
of mass~$c$; for the Einstein--Maxwell constraints, we match the
electric charge~$Q$. For large enough $R$ the patched data is
sufficiently close to flat, and we solve the constraints up to a
finite-dimensional obstruction space. To accomplish this, we iterate
linear corrections, in spaces with suitable decay at the boundary of
the annulus. The procedure is controlled using basic elliptic
estimates, both integral and pointwise, and with suitable weights. The
integral estimates for the time-symmetric case are relatively
straightforward, while the weighted Schauder estimates follow from the
standard interior estimates and scaling. From here, a topological
argument using the finite-dimensional space of parameters allows us to
cover the cokernel and solve the constraints. Similar comments apply
in the non-time-symmetric case.

The main first step, then, is how to solve the relevant linear
equations with compact support. For \xsmath{g= \r{g} + h} with $h$
small, the scalar curvature to leading order is governed by a
divergence:
{$$
R(\r{g} + h) = DR_{\r{g}} (h) + \mathcal{Q}(h)=
\mathrm{div}_{\r{g}}\mathrm{div}_{\r{g}}(h -
(\mathrm{tr}_{\r{g}} h) \r{g}) + \mathcal{Q}(h),
$$}\unskip
where $\mathcal{Q}(h)$ is a combination with bounded coefficients of
terms of components of $h\otimes \partial^2 h$ and $\partial h \otimes
\partial h$. A similar expansion holds for the constraints operator. At
the linear level, then, solving for a compactly supported perturbation
for the constraints map near \xsmath{(\r{g}, 0)} amounts to being able
to solve suitable divergence equations with compact support. 

Recent works~\cite{mao-tao, mao-oh-tao} have indeed developed this
approach near the flat data (with an extension to data near hyperbolic
data obtained in~\cite{ptc-cog-nut}). For instance, to solve a double
divergence equation, the authors use a Bogovski\v{i}-type operator as
follows. Let $\Omega\subset \mathbb{R}^3$ be an open set, which for
simplicity we take to be convex, and containing an open ball~$B$, and
let $\eta\in C^{\infty}_c(B)$ with $\int_{\mathbb{R}^3} \eta(x)\,
\mathrm{d}x=1$. Then for any $f\in C^\infty_c(\mathbb{R}^n)$ we define
the tensor $Sf$ to have Cartesian\break components
{$$
(Sf)^{ij}(x) = \int_{\mathbb{R}^3} \left(\int_{|x-y|}^\infty \eta
\left(r \frac{x-y}{|x-y|} + y\right) r^2\,\mathrm{d}r\right)
\frac{(x-y)^i(x-y)^j}{|x-y|^3}\, f(y)\,\mathrm{d}y.
$$}\unskip
If $f\in C^{\infty}_c(\Omega)$, then $Sf$ is supported in $\Omega$ as
well, and moreover, if $\int_{\Omega} f(x)\,\mathrm{d}x =0 =
\int_{\Omega} x^k f(x)\, \mathrm{d}x$, for $k\in \{1, 2, 3\}$, then
\xsmath{\mathrm{div}_{\r{g}} \mathrm{div}_{\r{g}} (Sf) = f}. 
The integral conditions are of course readily seen to be the
orthogonality of $f$ to the cokernel of the double divergence. One can
furthermore use a straightforward finite patching argument to handle
more general domains such as the annuli used in gluing constructions. 

As a remark related to this and to the Einstein--Maxwell equations, we
recall the well-known procedure for solving the divergence equation
\xsmath{\mathrm{div}_{\r{g}} X = f}, with $f$ and the solution $X$
compactly supported (see, e.g.\ \cite[Section~1.3]{nirenberg-fa}). Of
course we assume that $\int_{\mathbb{R}^n} f(x)\, \mathrm{d}x =0$. In
one dimension this is straightforward, as $X(x) = \int_{-\infty}^x
f(t)\, \mathrm{d}t$ satisfies the required conditions. We can proceed
inductively. Suppose $f$ is compactly supported in $D= D_n \times
I\subset \mathbb{R}^{n+1}$, where $D_n$ is an open rectangle in
$\mathbb{R}^n$ and $I\subset \mathbb{R}$ is an open interval. For $x\in
\mathbb{R}^n$, we let $F(x)= \int_{\mathbb{R}} f(x, t)\, \mathrm{d}t$.
By assumption, then $\int_{\mathbb{R}^n} F(x)\, \mathrm{d}x =0$, so by
induction, there are $Y^1, \ldots, Y^n$ supported in $D_n$ with
$\partial_i Y^i(x) = F(x)$. If~$\theta$ is supported in $I$ with
$\int_{\mathbb{R}} \theta (t)\, \mathrm{d}t =1$, we let $X^i(x,
x^{n+1}) = Y^i (x) \theta(x^{n+1})$ for\break $i\in \{1, \ldots, n\}$ and
{$$
X^{n+1}(x, x^{n+1}) = \int_{-\infty}^{x^{n+1}} (f(x, t) -
\theta (t) F(x))\,\mathrm{d}t.
$$}\unskip
The vector field $X$ with components $X^1, \ldots, X^{n+1}$ is
supported in $D$ and has divergence~$f$. 

Returning to the constraints, one would still need to solve the
nonlinear projected problem, which is essentially done by an iteration
scheme. That said, one major benefit of the Bogovski\v{i}-type
operator is that regularity is maintained when the data is in
appropriate Sobolev spaces, a question which arises when considering
the natural regularity for initial data for the evolution problem. 
Moreover the method has been extended beyond the asymptotically flat
regime to give an asymptotically hyperboloidal
construction~\cite{ptc-cog-nut}, and this is certainly intriguing since
the asymptotically hyperboloidal case proves to be more challenging
(for instance the conformal infinity is not a point as it is for the
asymptotically Euclidean case), see for example~\cite{ptc-delay-AH-cag,
ils-AH, ails-SF}. The Bogovski\v{i}-type operator method has also been
used for other localization results such as Carlotto--Schoen type
gravitational shielding constructions~\cite{mao-tao}. Finally, one
often wants to connect the solutions of the Einstein constraints to the
spacetimes they determine, whether in theory or numerically. 
Numerical implementation of gluing constructions seems quite
challenging, and though there have been a handful of works in this
direction~\cite{Giul-holz, Doulis-Rinne, dasz-frau, pook-giul},
alternative methods to do gluing constructions could prove useful if
they are more readily implemented numerically. It is natural to wonder
if the Bogovski\v{i}-type operator approach of~\cite{mao-oh-tao} might
be helpful in this direction, and likewise for the characteristic
gluing constructions of~\cite{aretakis-czimek-rodnianski} involving
transport-elliptic differential equations, which are of a different
character than the spacelike gluing.

\section{Further applications to spacetimes} \label{sec:further}

As mentioned earlier, one of the early goals for the gluing
construction was to create vacuum initial data which would evolve to a
vacuum spacetime that is asymptotically simply in the sense of Penrose
conformal compactification, a {\it purely radiative} spacetime. Since
work of Friedrich~\cite{friedrich} the following was known: if one
could construct a family of solutions to the vacuum constraints with
small mass tending to zero, each of which is Schwarzschild in a {\it
uniform} neighborhood of spacelike infinity, then for members of the
family with sufficiently small mass, the spacetime evolution would be
an asymptotically simple vacuum spacetime. The result
in~\cite{cor:schw} on its own did not quite give such a construction. 
The main issue, as one can see from (\ref{eq:mc}), is that as the mass
$m$ tends to zero, the projection onto the linear directions is on the
order of $mc$, and so for small mass, the center of mass may need to
grow larger to accommodate. To control the center, Chru\'sciel and
Delay~\cite{cd:as-pen} restrict to a family of parity-symmetric
metrics, for which the center of mass in the construction can always be
taken to vanish. In~\cite{cor:as-pen}, we construct a larger family of
solutions to the constraints with small masses and controlled centers,
by using the second variation of the mass at the Euclidean
data~\cite{cb-mar-pmt, brill-deser-pmt} to control the mass from below
in terms of the center. With either approach, the gluing to
Schwarzschild at infinity then finishes the initial data construction. 

In a completely different direction, Li and Yu~\cite{li-yu} construct
vacuum initial data, free of trapped surfaces, whose evolution will
contain a trapped surface (cf.~\cite{li-mei}). Of note, the solution
contains an interior region of Minkowskian data, while the exterior is
an exterior of a space-time slice of a Kerr spacetime. The construction
builds on the monumental work of Christodoulou on formation of trapped
surfaces {\it in vacuum}~\cite{chr-bh}, from data on null
hypersurfaces; in particular, the construction in~\cite{li-yu} uses
data given at past null infinity. 

While we do not have the space and time to discuss any details here,
there is one feature of the proof that we want to note. The
construction uses gluing to Kerr exterior as in~\cite{cs:ak}, but in a
regime where the limiting data is not Minkowksian, but Kerrian. As
such, the cokernel has {\it smaller} dimension than at the flat data,
so the gluing is not quite as obstructed. This has been observed in
other places, such as in initial data constructions with positive
cosmological constant~\cite{ptc-poll-kds} and in more recent work on
{\it multi-localized} initial data~\cite{acp:multi}. In the latter
work, the relevant approximate kernel is that of a Schwarzschild, which
is one-dimensional (for non-zero mass). The motivation for the
multi-localized initial data constructed in~\cite{acp:multi} comes from
the nonlinear stability result of Anderson and
Pasqualotto~\cite{AndPas22}. The template for the construction is the
Brill--Lindquist metric (\ref{eq:BLm}), and we glue ends on around each
of the punctures in the set $\mathfrak{p}$ to make initial data
on~$\mathbb{R}^3$. Of note is that we are able to handle certain
configurations where $\mathfrak{p}$ is infinite, and in particular, we
construct smooth metrics on $\mathbb{R}^3$ with zero scalar curvature
and infinitely many minimal spheres, each of them exterior to all the
others. 

As remarked earlier, in recent years {\it characteristic gluing} has
been developed in a series of papers by Aretakis, Czimek and
Rodnianski; for an overview and further references,
see~\cite{aretakis-czimek-rodnianski}. For the characteristic initial
value problem, one specifies certain free data on what will be two
transversely intersecting null hypersurfaces in the vacuum spacetime,
after which one can appeal to local well-posedness results for the
characteristic initial value problem for the Einstein vacuum
equation~\cite{rend-char, luk-char}. For the characteristic gluing
problem, one seeks to glue data on two spheres along a null
hypersurface. As with the spacelike case, the relevant equations that
must be satisfied arise from the vacuum Einstein equation together with
the fundamental compatibility equations for an embedding, in this case
the compatibility coming from a double null foliation. The resulting
equations for metric and curvature components are the {\it null
structure equations}, of transport-elliptic character along either the
ingoing or outgoing null hypersurface. As is noted
in~\cite{aretakis-czimek-rodnianski}, there are obstructions to gluing
data from two spheres, such as a monotonicity deriving from the
Raychaudhuri equation.

There are a number of noteworthy applications of their technique, and
we just mention a few. Their methods can recover the asymptotic gluing
to Kerr as in~\cite{cs:ak}, as well as localized gluing
from~\cite{car-sch}, and without loss of decay in the transition
region. In~\cite{Czi-Rod}, the method is employed to achieve {\it
obstruction-free} gluing. The asymptotic gluing construction along the
lines of~\cite{cor:schw, cs:ak} glues initial data with certain
energy--momenta to Kerr data with energy--momenta which is close to the
starting one, as we say above, tuned to cover the linear obstruction
from the approximate cokernel. A natural question is whether one can
glue pieces from spacetimes whose parameters are not necessarily close,
and results in this direction have been obtained in~\cite{Czi-Rod}. 

The characteristic gluing construction has been employed by Kehle and
Unger~\cite{Kehle-Unger} to perform a specific {\it event horizon
gluing} in the context of the Einstein--Maxwell charged scalar field
system, in spherical symmetry. As an upshot, they produce a
counterexample to the third law of black hole thermodynamics. 

\section{Conclusion} \label{sec5}

Gluing methods for constructing solutions of the Einstein constraint
equations have garnered a number of noteworthy achievements. While
more than twenty years have passed since the gluing results
of~\cite{cor:schw} and~\cite{imp}, new techniques and applications are
still being developed and continue to highlight the interplay between
initial data results and the evolution problem. 

\section*{Dedication}

This article is dedicated to Mme. Yvonne Choquet-Bruhat, for giving so
much to so many of us for so many years. The memory of Yvonne that is
emblazoned upon me comes from over twenty years ago at a workshop in
Carg\`ese, Corsica, marking the fiftieth anniversary of her landmark
paper on the initial-value problem. Walking back to the seminar room
with a colleague, we encountered Yvonne taking a bit of respite from
the sunshine under a tree; while it was not an apple tree, what sprung
to my imagination was the universe (or the spirit of Sir Isaac Newton)
willing one to fall and land in the palm of her hand.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\section*{Acknowledgments}

The author would like to thank Farhan Abedin, David A.~Maxwell and Iva
Stavrov Allen for their comments on a draft manuscript. The author is
supported in part by the AMS-Simons Research Enhancement Grant for PUI
Faculty.

\back{}

\bibliographystyle{crunsrt}
\bibliography{crmeca20240924}
\refinput{crmeca20240924-reference.tex}

\end{document}
