\section{Introduction}

One of the most striking evidence for the formation of a quark-gluon plasma in heavy-ion collisions at RHIC and at the LHC is the suppression of high $p_T$ jets, phenomenon known as jet quenching. If the measurement of such suppression is of fundamental importance \textit{per se}, it is not enough to constrain the various theoretical models for plasma-jet interactions. It is then crucial to confront models with other jet observables. Among them, jet substructure observables are promising probes of the modifications of the partonic cascades produced in the medium \cite{Andrews:2018jcm}.
In this paper, we discuss one such observable: the jet fragmentation function.

The theoretical framework for this study is perturbative QCD (pQCD) and its emerging picture for jet evolution detailed in \cite{Caucal:2018dla,Caucal:2019uvr}. This picture includes two kinds of radiation: standard vacuum-like emissions (VLEs) triggered by the parton virtuality and medium-induced emissions (MIEs) triggered by multiple collisions inside the medium, characterised by its quenching parameter $\hat{q}$.
We have shown in \cite{Caucal:2018dla} that the partonic cascades can be factorised in three steps: {\tt (i)} a cascade of VLEs {\em inside the
    medium} satisfying $\omega^3\theta^4>2\hat{q}$ and $\theta>\theta_c \equiv
  2/\sqrt{{\hat q}L^3}$ with $\omega$ and $\theta$ the energy and angle of emission and $L$ the jet path length,
  %
{\tt (ii)} each parton escaping the previous vacuum-like cascade propagates
  through the medium over a distance $L$ and thus sources
  MIEs,
  %
{\tt (iii)} finally all the outgoing partons trigger another cascade of VLEs outside
  the medium in the phase space $t_f=2/\omega\theta^2>L$. The first emission outside the medium is not constrained by angular ordering because of the colour decoherence due to the interactions with the plasma \cite{MehtarTani:2010ma,MehtarTani:2011tz,CasalderreySolana:2011rz}.
The implementation of this picture within a Monte-Carlo (MC) parton shower is straightforward (see \cite{Caucal:2018ofz,Caucal:2019uvr} for detailed descriptions). All the numerical results presented here are obtained using this~MC.

This paper is organised as follows: we first present our results for the ``standard'' fragmentation function defined from the hadrons inside jets. This observable being infrared and collinear (IRC) unsafe, we define a new fragmentation function under better theoretical control in pQCD. This IRC safe fragmentation function is then studied in the vacuum and in nucleus-nucleus collisions.

\section{Nuclear modifications of the standard fragmentation function}

The standard jet fragmentation function $\mathcal{D}(x)$ and its nuclear modification factor $\mathcal{R}(x)$ are defined as
\begin{equation}\label{frag-def}
 \mathcal{D}(x)=\frac{1}{N_{\textrm{jets}}}\frac{\textrm{d} N}{\textrm{d} x},\qquad \mathcal{R}(x)=\frac{\mathcal{D}^{\textrm{med}}(x)}{\mathcal{D}^{\textrm{vac}}(x)}
 \end{equation}
where $N_{\textrm{jets}}$ is the number of selected jets and $\textrm{d} N$ is the number of jet constituents (hadrons) with a momentum fraction $x$ between $x$ and $x+\textrm{d} x$. The integral of $\mathcal{D}(x)$ over $x$ is equal to the mean total intrajet multiplicity of hadrons, a quantity which is of course not IRC safe. Thus, the fragmentation function is strongly sensitive to hadronisation corrections
or on the unphysical cut-off $k_{\perp,\rm min}$ required to regulate the collinear divergence in a parton level calculation.
\begin{figure}[t]
  \centering
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[page=1,width=\textwidth]{ff-nonphys.pdf}
    \caption{\small Variations in  $\theta_{\rm max}$ and $k_{\perp,\text{min}}$.}\label{Fig:MCunphys}
  \end{subfigure}
  %
  \hfill
  %
  \begin{subfigure}[t]{0.48\textwidth}
    \includegraphics[page=1,width=\textwidth]{ff-phys.pdf}
    \caption{\small Variations in  $\hat{q}$, $L$ and $\alpha_{s,\rm med}$.}
    \label{Fig:MC-pheno}
  \end{subfigure}
  \caption{\small (Left) Variability of our MC results for the ratio $\mathcal{R}(x)$ w.r.t. changes
  in $k_{\perp,\rm min}$ and the maximal branching angle $\theta_{\rm max}$. The envelop is dominated by $k_{\perp,\rm min}$ variations. (Right) MC results for the 4 sets of values of
  the medium parameters $\hat{q}$, $L$ and the strong coupling for medium-induced vertices $\alpha_{s,\rm med}$ that provide similarly good descriptions
  of the LHC data~\cite{Aaboud:2018twu} for the nuclear  modification factor for jets $R_{AA}$~\cite{Caucal:2019uvr}.}
  \label{fig-R-param-dependence}
\end{figure}

To illustrate this, we show on Fig.~\ref{fig-R-param-dependence}-left the variability of $\mathcal{R}(x)$ with respect to the cut-off $k_{\perp,\rm min}$ as given by our MC, which does not currently include any hadronisation model. As expected, the envelop of the variations is large. Nevertheless, the general trend of the curve is in qualitative agreement with the recent ATLAS data \cite{Aaboud:2018hpb}. Namely, one observes a significant enhancement of the ratio at small $x$ and for $x$ close to $1$, while for $0.02\lesssim x \lesssim 0.2$, the ratio becomes smaller than~$1$. We point out that since both $x\mathcal{D}^{\textrm{med}}(x)$ and $x\mathcal{D}^{\textrm{vac}}(x)$ are normalized to~1, one cannot have $\mathcal{R}(x)>1$ for all $x$. That said, we now briefly explain the physical mechanisms leading to these enhancements at both ends of the spectrum.

In the small $x$ part, these mechanisms differ significantly from those encountered in the literature~\cite{Casalderrey-Solana:2016jvj,Tachibana:2017syd,KunnawalkamElayavalli:2017hxo},
where the increase of soft particle production at
relatively large angles from the jet axis is related to the medium response to the jet propagation. In our picture, soft particle production is dominated by gluons produced outside the medium with $t_f>L$ at the end of the evolution. This evolution is amplified in the presence of a medium w.r.t.\ the vacuum because of \cite{Caucal:2020xad}: (i) the violation of angular ordering by the first emission outside the medium, which
re-opens the angular phase-space for subsequent evolution \cite{Mehtar-Tani:2014yea,Caucal:2018dla}, and (ii) the presence
of relativitely hard MIEs with $\omega\gg\omega_{\rm br}\equiv \alpha_{s}^2\hat{q}L^2$ which remain inside the jet and can further radiate VLEs outside the medium~\cite{Caucal:2019uvr}. (The multiple branching scale $\omega_{\rm br}$ is the typical scale for energy loss at large angles, i.e.\ MIEs deviated out of the jet cone \cite{Blaizot:2015lma}.)

For $x$ close to 1, we have shown in \cite{Caucal:2020xad} that the enhancement seen in the ratio is essentially \textit{not} a change in the in-medium fragmentation process but rather a bias toward hard fragmenting jets induced by the steeply falling jet spectrum.
This spectrum favors jets that lose less energy. In hard-fragmenting jets, the leading parton\footnote{Most likely a quark, since quark jets have a harder core and lose less energy than gluon jets \cite{Spousta:2015fca,Caucal:2020xad}.} sourcing the jet must not have radiated too many partons in order to have a final $x$ fraction close to 1. As the intrajet multiplicity is positively correlated with the energy loss \cite{Caucal:2019uvr}, hard fragmenting jets lose less energy than average jets, hence are favored in nucleus-nucleus collisions.
This argument is actually very general as it applies to many theoretical models \cite{KunnawalkamElayavalli:2017hxo,Rajagopal:2016uip,Casalderrey-Solana:2016jvj,Casalderrey-Solana:2018wrw}. It also explains the strong correlation between the nuclear modification factor for jets $R_{AA}$ and the $x\sim 1$ region of the fragmentation function, since both quantities are controlled by energy loss effects. This correlation is apparent on Fig.~\ref{fig-R-param-dependence}-right: the 4 sets of physical parameters in our calculations that all give a good description of $R_{AA}$ \cite{Caucal:2019uvr} scale identically for $x\sim1$ whereas the degeneracy is lifted for $x\ll 1$. We conclude that by studying the soft sector in the fragmentation of jets, one can better constrain theoretical models that have been tuned to yield a good description of $R_{AA}$.

\section{Fragmentation into subjets}

In this section, we first define an IRC safe fragmentation function relying on subjets and we calculate it at leading logarithmic (LL) accuracy in the vacuum. Then, we argue that this new observable is sensitive to the same medium effects while being resilient to non-perturbative physics.

\paragraph{Definition and LL result in the vacuum.}The basic idea of this observable is to replace the unphysical cut-off $k_{\perp,\rm min}$ by a parameter of the observable, $k_{\perp,\rm cut}$ which is chosen much larger than the confinement scale $\Lambda_{\rm QCD}$. For a given jet with transverse momentum $p_{T,\textrm{jet}}$, we decluster the jet using the Cambridge/Aachen algorithm \cite{Dokshitzer:1997in,Wobisch:1998wt} following the hardest branch. This gives a list of subjets. In this list, we keep only those which have a transverse momentum with respect to the hard branch larger than $k_{\perp,\textrm{cut}}$. The subjet fragmentation  is then defined as the number $\textrm{d} N_{\textrm{sub}}$ of subjets with \textit{splitting fraction} $z<1/2$ between $z$ and $z+\textrm{d} z$ normalised by the total number of selected jets~\cite{Caucal:2020xad}:
\begin{equation}
 \mathcal{D}_{\textrm{sub}}(z)\equiv\frac{1}{N_\textrm{jets}}\frac{\textrm{d} N_{\textrm{sub}}}{\textrm{d} z}
\end{equation}
This definition is close to the Iterated Soft Drop multiplicity differential in $z$ \cite{Frye:2017yrw}.

To calculate this quantity in the vacuum, it is convenient to introduce the generating functional of the exclusive probability distribution $P _i^{(n)}(z_1,...,z_n)$ of $n$ subjets with splitting fractions $z_1$,...,$z_n$ in a jet with transverse momentum $p_T$, opening angle $\theta$ and ``flavor'' $i\in\{q,g\}$:
 \begin{equation}\label{Zjet}
 Z_i(p_T,\theta|u(z))=\sum_{n=1}^{\infty}\int \textrm{d} z_1...\textrm{d} z_n\, u(z_1)...u(z_n)P_i^{(n)}(z_1,...,z_n)
\end{equation}
The fragmentation function into subjets $\mathcal{D}_{i,\rm sub}(z|p_T,R)$ of a $i$-jet with opening angle $R$ and transverse momentum $p_T$ is then obtained by a functional differentiation of $Z_i(p_T,R|u(z))$ with respect to the probing function $u(z)$, at $u=1$.
From the Markovian property of the branching process ordered in angles,  $Z_i$ satisfies the following differential equation in the vacuum:
\begin{equation}\label{Z-primary}
\frac{\partial Z^{\rm vac}_{i}(p_T,\theta|u)}{\partial \log(\theta)}=\sum_{(a,b)}\int_0^{1/2}\textrm{d} z\,\frac{\alpha_s(k_\perp)}{\pi}\Theta_{\rm cut}\Phi_i^{ab}(z)\left(u(z)Z^{\rm vac}_b((1-z)p_T,\theta|u)-Z^{\rm vac}_i(p_T,\theta|u)\right)
\end{equation}
where $\Phi_i^{ab}(z)$ are the unregularised DGLAP splitting functions (the sum runs over all distincts pairs of partons $(a,b)$) and $\Theta_{\rm cut}$ is a step function enforcing $k_\perp>k_{\perp,\rm cut}$. Eq.~\ref{Z-primary} includes corrections associated with the recoil of the hard branch. At LL accuracy, one can ignore such contributions by setting $Z_b((1-z)p_T,\theta|u)\simeq Z_b(p_T,\theta|u)$ and one can approximate $\Phi_i^{ab}(z)\simeq 2C_i/z$ and $k_\perp\simeq zp_T\theta$ in the argument of $\alpha_s$. At this stage, it is possible to solve exactly the corresponding evolution equation for $\mathcal{D}_{\rm sub}(z|p_T,\theta)$ with initial condition $\mathcal{D}_{\rm sub}(z|p_T,0)=0$. Furthermore, if one neglects quark-gluon mixing\footnote{Even though it formally matters at LL accuracy.}, one finds a good approximation of the LL result under the simple form:
\begin{equation}\label{frag-sub-vac}
\mathcal{D}^{\rm vac}_{i,\rm sub}(z|p_T,R)\simeq \frac{2C_i}
{\pi}\int_{0}^R\frac{\textrm{d}\theta}{\theta}\,
\frac{\alpha_s(zp_T\theta)}{z}\Theta(zp_T\theta-k_{\perp,\rm cut})
\end{equation}
This formula shows explicitly how the parameter $k_{\perp,\rm cut}\gg \Lambda_{\rm QCD}$ plays also the role of an infrared regulator in the $\theta$ integral. It makes manifest the connection between $\mathcal{D}_{\rm sub}(z)$ and the primary Lund plane density integrated over all angles with $z\simeq k_\perp/(\theta p_{T,\textrm{jet}})$ fixed \cite{Dreyer:2018nbf,Lifson:2020gua}.  Finally, the resummed result \eqref{frag-sub-vac} should be matched with a fixed order calculation for quantitative comparisons with data in $pp$ collisions.

\begin{figure}[t]
  \centering
  \includegraphics[page=2,width=0.48\textwidth]{plot-subjets.pdf}
  \hfill
   \includegraphics[page=5,width=0.48\textwidth]{plot-subjets.pdf}
  \caption{\small Variability of our MC results for the ratio $\mathcal{R}_{\rm sub}(z)$ w.r.t.\ changes
  in $k_{\perp,\rm min}$ and  $\theta_{\rm max}$ (left) and w.r.t\ the medium parameters (right) that reproduce the ATLAS $R_{AA}$ ratio as in Fig.~\ref{Fig:MC-pheno}. For the latter figure, $k_{\perp,\textrm{cut}}=2$~GeV and the unphysical parameters are fixed to
$\theta_{\textrm{max}}=1$ and $k_{\perp,\rm min}=250$~MeV.
   }
\label{Fig:subjet-med}
\end{figure}

\paragraph{Monte-Carlo results in Pb+Pb collisions.}We now discuss the nuclear modification factor $\mathcal{R}_{\rm sub}(z)$ defined as in \eqref{frag-def} and calculated with our MC. First of all, being IRC safe, $\mathcal{R}_{\rm sub}(z)$ is much less sensitive to $k_{\perp,\rm min}$ as shown Fig.~\ref{Fig:subjet-med}-left. The larger $k_{\perp,\rm cut}$ is, the more this sensitivity is reduced. On Fig.~\ref{Fig:subjet-med}-right, the ratio $\mathcal{R}_{\rm sub}(z)$ is shown for the 4 sets of parameters that fit $R_{AA}$ measured by ATLAS~\cite{Aaboud:2018twu}. The reduction at large $z\sim0.5$ is again essentially a consequence of the normalisation factor $N_{\rm jets}$ \cite{Caucal:2020xad}. Jets with $z\sim 0.5$ and $k_\perp>k_{\perp,\rm cut}$ are typically two-prongs jets with angular separation larger than the coherence angle $\theta_c\propto(\hat{q}L^3)^{-1/2}$, and thus lose \textit{more} energy than average jets \cite{Caucal:2019uvr,Casalderrey-Solana:2019ubu}. They are disfavored by the initial jet spectrum and suppressed in Pb-Pb collisions. On the other hand, the increasing behaviour of the ratio as $z$ decreases is mainly a consequence of intrajet semi-hard MIEs \cite{Caucal:2020xad}. Such emissions source subjets that are finally captured in the declustering procedure performed to calculate $\mathcal{D}_{\textrm{sub}}^{\rm med}(z)$. As $k_{\perp,\rm cut}$ is much larger than the scale $2/(LR)$, $\mathcal{D}_{\textrm{sub}}^{\rm med}(z)$ probes only a tiny portion of the ``outside'' region of phase space with $t_f>L$. Angular ordering violation has therefore a less pronounced effect than for the standard ratio $\mathcal{R}(x)$.

To conclude, we point out that $k_{\perp,\rm cut}$ has a determinant role in heavy-ion collisions. It must not be too small to avoid large hadronisation corrections. However, if it is too large, the interesting medium effects in the relatively soft sector highlighted above are not captured. In that respect, substructure observables relying on dynamical grooming is another interesting possibility \cite{Mehtar-Tani:2019rrk,Soto-Ontoso:2020ola}.

\paragraph{Acknowledgements} The work of P.C.,  E.I. and G.S. is supported in part by the Agence Nationale de la Recherche project
 ANR-16-CE31-0019-01.   The work of A.H.M.
is supported in part by the U.S. Department of Energy Grant \# DE-FG02-92ER40699.
\bibliographystyle{JHEP}


\subsection{Laplace obstacle problem}

Let $\Omega \subset \mathbb{R}^d$, $d \ge 2$, be a bounded Lipschitz domain. Let $\chi$ be an affine function
with $\chi \le 0$ on $\partial\Omega$. Denote the set of admissible functions by
\begin{align}
 \mathcal{A} := \set{v \in \mathcal{X}}{v \ge \chi \text{ a.e.\ on } \Omega} \quad\text{with }\mathcal{X}:=H^1_0(\Omega)
\end{align}
and note that $\mathcal{A} \neq \emptyset$.
The minimization problem reads: Given a continuous linear functional $f \in H^{-1}(\Omega)$ with , find $u \in \mathcal{A}$ such that
\begin{align}\label{eq:obstacle:strong}
 E(u) = \min_{v \in \mathcal{A}} E(v),
 \quad \text{where }
 E(v) := \frac{1}{2} \norm{\nabla v}{L^2(\Omega)}^2 - \int_\Omega f v \d{x}.
\end{align}
Since $\mathcal{A} \neq \emptyset$ is convex and closed, it is well-known~\cite[Theorem~II.2.1]{MR1786735} that the minimization problem~\eqref{eq:obstacle:strong} admits a unique solution $u \in \mathcal{A}$.

We consider regular triangulations $\mathcal{T}_H$ of $\Omega$ into non-degenerate compact simplices and the corresponding first-order Courant finite element space
\begin{align}\label{eq:courant}
 \mathcal{X}_H := \set{v_H \in H^1_0(\Omega)}{v_H|_T \text{ is affine for all } T \in \mathcal{T}_H}.
\end{align}
Then, $\mathcal{A}_H := \mathcal{A} \cap \mathcal{X}_H \neq \emptyset$ is closed and convex. As in the continuous case, there hence exists a unique discrete minimizer $u_H \in \mathcal{A}_H$ such that
\begin{align}
 E(u_H) = \min_{v_H \in \mathcal{A}_H} E(v_H).
\end{align}
Under additional regularity $f \in L^2(\Omega)$, one can argue as in~\cite{bch2007} to show reliability
\begin{align}\label{eq:obstacle:reliable}
 \frac{1}{2} \norm{\nabla(u - u_H)}{L^2(\Omega)}^2
 \le E(u_H) - E(u)
 \le C_{\conststyle rel}^2 \eta_H^2,
\end{align}
for the residual error estimator with local contributions
\begin{align}\label{eq:obstacle:estimator}
 \eta_H(T)^2 = h_T \norm{\llbracket \partial_n u_H \rrbracket}{L^2(\partial T\cap\Omega)}^2
 + h_T^2 \sum_{\substack{E \in \mathcal{E}_H(T) \\ E \subset \partial\Omega}} \norm{f}{L^2(T)}^2
 + h_T^2 \sum_{\substack{E \in \mathcal{E}_H(T) \\ E \not\subset\partial\Omega}} \norm{f - f_E}{L^2(T)}^2,
\end{align}
where $h_T := |T|^{1/d}$ and $\mathcal{E}_H(T)$ is the set of all $(d-1)$-dimensional facets (i.e., edges for $d = 2$) and $f_E \in \mathbb{R}$ is the integral mean of $f$ over the corresponding patch $\omega_E = T \cup T'$ with $E = T \cap T'$. In addition to~\eqref{eq:union_of_sons}, we suppose uniform contraction of the mesh-size on refined elements, i.e., there exists $0 < \contract{ctr} < 1$
\begin{align}\label{eq:reduction:h}
 |T'| \le \contract{ctr} |T|
 \quad \text{for all \ } T \in \mathcal{T}_H \in \mathbb{T} \text{ \ and all \ } T' \in \mathcal{T}_h \in \mathbb{T}(\mathcal{T}_H) \text{ \ with \ } T' \subsetneqq T.
\end{align}
Then, stability~\eqref{eq:axiom:stability} and reduction~\eqref{eq:axiom:reduction} with $S(t) = C_{\conststyle stab} \,t$ and $R(t) = C_{\conststyle red} \, t^2$ follow as for the linear case; see~\cite{ckns2008} (or~\cite{pp2011} for the obstacle problem). The constants $C_{\conststyle rel}, C_{\conststyle stab} > 0$ depend only on $\Omega$, $d$, and uniform $\gamma$-shape regularity of the admissible meshes $\mathcal{T}_H \in \mathbb{T}$ in the sense of
\begin{align}\label{eq:gamma-shape-regular}
 \gamma := \sup_{\mathcal{T}_H \in \mathbb{T}} \max_{T \in \mathcal{T}_H} \frac{{\rm diam}(T)}{|T|^{1/d}} < \infty,
\end{align}
while $C_{\conststyle red} > 0$ and $0 < \contract{red} < 1$ depend additionally on $\contract{ctr}$.
The {\sl a~priori} convergence~\eqref{assumption:limit_u} follows essentially as in the seminal work~\cite{bv1984}: The assumption~\eqref{eq:union_of_sons} on the mesh-refinement implies that refinement leads to nested spaces, i.e., Algorithm~\ref{alg:abstract_algorithm} leads to $\mathcal{X}_\ell \subseteq \mathcal{X}_{\ell+1}$ and hence $\mathcal{A}_\ell \subseteq \mathcal{A}_{\ell+1}$ for all $\ell \in \mathbb{N}_0$. Therefore, $\mathcal{A}_\infty := {\rm closure}\big(\bigcup_{\ell \in \mathbb{N}_0} \mathcal{A}_\ell\big) \neq \emptyset$ is a closed and convex subset of $\mathcal{X}$ and thus gives rise to a unique minimizer $u_\infty \in \mathcal{A}_\infty$ such that
\begin{align}
 E(u_\infty) = \min_{v_\infty \in \mathcal{A}_\infty} E(v_\infty).
\end{align}
Based on estimates for the equivalent variational inequalities (see~\cite[Theorem~II.2.1]{MR1786735}), it follows that
\begin{align}
 \norm{u_\infty - u_\ell}{\mathcal{X}}^2
 \lesssim \inf_{v_\ell \in \mathcal{A}_\ell}
 \norm{u_\infty - v_\ell}{\mathcal{X}}
 \to 0 \quad \text{as } \ell \to \infty,
\end{align}
where we stress the different powers of the norms which are due to lack of Galerkin orthogonality; see, e.g.,~\cite{falk1974}. Overall, we thus get the following plain convergence result, where we note that for D\"orfler marking~\eqref{eq:marking:doerfler} with sufficiently small marking parameter $0 < \theta \ll 1$, \cite{ch2015} proves even rate-optimal convergence for $d = 2$.

\begin{proposition}
As long as the mesh-refinement strategy guarantees regular simplicial triangulations satisfying~\eqref{eq:union_of_sons}, \eqref{eq:reduction:h}, and~\eqref{eq:gamma-shape-regular}
 and as long as the marking strategy satisfies~\eqref{eq:marking}, Algorithm~\ref{alg:abstract_algorithm} for the Laplace obstacle problem~\eqref{eq:obstacle:strong} driven by the indicators~\eqref{eq:obstacle:estimator} yields convergence
\begin{align}\tag*{\qed}
 C_{\conststyle rel}^{-1} \, \norm{\nabla(u - u_\ell)}{L^2(\Omega)}
 \le \eta_\ell \to 0
 \quad \text{as } \ell \to \infty.
\end{align}
\end{proposition}

\subsection{Fractional Laplacian}\label{section:fractional}

Let $\Omega \subset \mathbb{R}^d$ be a bounded Lipschitz domain, $d \ge 2$, and $0 < s < 1$. Given $f \in H^{-s}(\Omega)$, we consider the Dirichlet problem of the fractional Laplacian
\begin{align}\label{eq:fractional:strong}
 (-\Delta)^s u = f \text{ in } \Omega
 \quad\text{subject to}\quad
  u = 0 \text{ in } \mathbb{R}^d \backslash \overline\Omega.
\end{align}
There are several different ways to define
$(-\Delta)^s$, e.g., in terms of the Fourier transformation~\cite{MR3893441}, via semi-group theory~\cite{MR3613319}, or as Dirichlet-to-Neumann map of a half-space extension problem~\cite{MR2354493}. For the latter, a convenient representation of the fractional Laplacian is given in terms of a principal value integral
\begin{align}
\big( (-\Delta)^s u \big)(x) := C(d,s) \,\, {\rm p.v.}\!\!\! \int_{\mathbb{R}^d} \frac{u(x) - u(y)}{|x-y|^{d+2s}} \d{y}
 \text{ \ with \ }
 C(d,s) := -2^{2s} \, \frac{\Gamma(s+d/2)}{\pi^{d/2}\Gamma(-s)},
\end{align}
where $\Gamma(\cdot)$ denotes the Gamma function. According to~\cite[Theorem~1.1]{MR3613319}, the weak formulation of~\eqref{eq:fractional:strong} reads: Find $u \in \mathcal{X} := \widetilde H^s(\Omega)$ such that
\begin{align}\label{eq:fractional:weak}
 a(u,v) := \frac{C(d,s)}{2} \int\!\!\!\!\int_{\mathbb{R}^d\times\mathbb{R}^d} \!\!\!\!\!\!\!\!\! \frac{[u(x)-u(y)][v(x)-v(y)]}{|x-y|^{d+2s}} \d{x}\d{y}
 = \int_\Omega fv \d{x}
 \quad \text{for all } v \in \mathcal{X}.
\end{align}
The Lax--Milgram lemma proves existence and uniqueness of $u \in \mathcal{X}$.

Following~\cite{MR3725761,fmp2019+}, we consider regular triangulations $\mathcal{T}_H$ of $\Omega$ into non-degenerate compact simplices and the corresponding first-order Courant finite element space $\mathcal{X}_H$ from~\eqref{eq:courant}.
Let $u_H \in \mathcal{X}_H$ the corresponding Galerkin solution, i.e.,
\begin{align}
 a(u_H, v_H) = \int_\Omega f v_H \d{x}
 \quad \text{for all } v_H \in \mathcal{X}_H.
\end{align}
Under the additional regularity assumption $f \in L^2(\Omega)$, the local contributions of the error estimator from~\cite{fmp2019+} read
\begin{subequations}\label{eq:fractional:estimator}
\begin{align}
 \eta_H(T) :=  \, \norm{h_H^s \,[ f - (-\Delta)^s u_H] }{L^2(T)}
 \quad \text{for all } T \in \mathcal{T}_H,
\end{align}
with the modified local mesh-width
\begin{align}
 h_H^s|_T := \begin{cases}
 |T|^{s/2} & \text{for } 0 < s \le 1/2, \\
 |T|^{1/4} \dist(\cdot,\partial T)^{s-1/2} & \text{for } 1/2 \le s < 1.
 \end{cases}
\end{align}
\end{subequations}
According to~\cite[Theorem~2.3]{fmp2019+}, the error estimator is reliable
\begin{align}\label{eq:reliable}
 \norm{u - u_H}{\mathcal{X}} \le C_{\conststyle rel} \, \eta_H.
\end{align}
Provided~\eqref{eq:union_of_sons} and~\eqref{eq:reduction:h}, stability~\eqref{eq:axiom:stability} and reduction~\eqref{eq:axiom:reduction} are proved in~\cite[Proposition~3.1]{fmp2019+} with $S(t) = C_{\conststyle stab} \,t$ and $R(t) = C_{\conststyle red} \, t^2$. The constants $C_{\conststyle rel}, C_{\conststyle stab} > 0$ depend only on $\Omega$, $d$, $s$, and uniform $\gamma$-shape regularity~\eqref{eq:gamma-shape-regular}, while $C_{\conststyle red} > 0$ and $0 < \contract{red} < 1$ depend additionally on $\contract{ctr}$.
 Finally, {\sl a~priori} convergence~\eqref{assumption:limit_u} follows as in the seminal work~\cite{bv1984} (and essentially with the same arguments as in the previous section): The assumption~\eqref{eq:union_of_sons} on the mesh-refinement implies that refinement leads to nested spaces, i.e., Algorithm~\ref{alg:abstract_algorithm} leads to $\mathcal{X}_\ell \subseteq \mathcal{X}_{\ell+1}$ for all $\ell \in \mathbb{N}_0$. Therefore, $\mathcal{X}_\infty := {\rm closure}\big(\bigcup_{\ell \in \mathbb{N}_0} \mathcal{X}_\ell\big)$ is a closed subspace of $\mathcal{X}$ and the Lax--Milgram lemma guarantees existence and uniqueness of $u_\infty \in \mathcal{X}_\infty$ such that
\begin{align}
 a(u_\infty, v_\infty) = \int_\Omega fv_\infty \d{x}
 \quad \text{for all } v_\infty \in \mathcal{X}_\infty.
\end{align}
With the Galerkin orthogonality and the resulting C\'ea lemma, it follows that
\begin{align}
 \norm{u_\infty - u_\ell}{\mathcal{X}} \lesssim \min_{v_\ell \in \mathcal{X}_\ell} \norm{u_\infty - v_\ell}{\mathcal{X}}
 \to 0 \quad \text{as } \ell \to \infty.
\end{align}
Overall, we thus get the following plain convergence result, where we note that for D\"orfler marking~\eqref{eq:marking:doerfler} with sufficiently small marking parameter $0 < \theta \ll 1$, \cite[Theorem~2.6]{fmp2019+} proves even rate-optimal convergence.

\begin{proposition}
As long as the mesh-refinement strategy guarantees regular simplicial triangulations satisfying~\eqref{eq:union_of_sons}, \eqref{eq:reduction:h}, and~\eqref{eq:gamma-shape-regular}
 and as long as the marking strategy satisfies~\eqref{eq:marking}, Algorithm~\ref{alg:abstract_algorithm} for the fractional Laplacian~\eqref{eq:fractional:weak} driven by the indicators~\eqref{eq:fractional:estimator} yields convergence
\begin{align}\tag*{\qed}
 C_{\conststyle rel}^{-1} \, \norm{u - u_\ell}{\widetilde{H}^s(\Omega)}
 \le \eta_\ell \to 0
 \quad \text{as } \ell \to \infty.
\end{align}
\end{proposition}

\subsection{Weakly-singular integral equations}

Let $\Omega \subset \mathbb{R}^d$ with $d = 2,3$ be a Lipschitz domain with compact boundary $\partial\Omega$ such that $\diam(\Omega)<1$ if $d=2$.
On a (relatively) open subset $\Gamma \subseteq \partial\Omega$ and given $f \in H^{1/2}(\Gamma)$, the weakly-singular integral equation
\begin{align}\label{eq:weaksing:strong}
(Vu)(x) := \int_\Gamma G(x-y) u(y) \d{y} = f(x)
 \quad \text{for all } x \in \Gamma
\end{align}
seeks the unknown integral density $u \in \mathcal{X} := \widetilde{H}^{-1/2}(\Omega)$. Here, $G(\cdot)$ is the fundamental solution of the Laplacian, i.e., $G(z) = -\frac{1}{2\pi} \log|z|$ for $d = 2$ resp.\ $G(z) = \frac{1}{4\pi} |z|^{-1}$ for $d = 3$.
We note that for $\Gamma = \partial\Omega$, \eqref{eq:weaksing:strong} is equivalent to the Dirichlet problem
\begin{align*}
 -\Delta U = 0 \text{ in } \Omega
 \quad \text{subject to } U = f \text{ on } \partial\Omega,
\end{align*}
supplemented by the appropriate radiation condition if $\Omega$ is unbounded; see~\cite{mclean}.
The weak formulation reads
\begin{align}\label{eq:weaksing:weak}
 a(u, v) := \int\!\!\!\!\int_{\Gamma\times\Gamma} G(x-y) u(y) v(x) \d{y} \d{x}
 = \int_\Gamma fv \d{x}
 \quad \text{for all } v \in \mathcal{X},
\end{align}
and the Lax--Milgram lemma yields existence und uniqueness of the solution $u \in \mathcal{X}$.

For a fixed polynomial degree $p \ge 0$ and a regular triangulation $\mathcal{T}_H$ of $\Gamma$ into non-degenerate compact surface simplices, we consider standard boundary element spaces $\mathcal{X}_H = \mathcal P^p(\mathcal{T}_H)$ consisting of $\mathcal{T}_H$-piecewise polynomials of degree $\le p$ (w.r.t.\ the boundary parametrization). Let $u_H \in \mathcal{X}_H$ be the corresponding Galerkin solution, i.e.,
\begin{align}
 a(u_H, v_H) = \int_\Gamma f v_H \d{x}
 \quad \text{for all } v_H \in \mathcal{X}_H.
\end{align}
According to the seminal work~\cite{cms2001} and under additional regularity $f \in H^1(\Gamma)$, the local contributions of the residual error estimator read
\begin{align}\label{eq:weaksing:estimator}
 \eta_H(T) := h_T^{1/2} \, \norm{\nabla_\Gamma(f - Vu_H)}{L^2(T)}
 \quad\text{with}\quad
 h_T := |T|^{1/(d-1)},
\end{align}
where $\nabla_\Gamma(\cdot)$ is the surface gradient and $|\cdot|$ is the surface measure. While~\cite{cms2001} proves reliability~\eqref{eq:reliable}, stability~\eqref{eq:axiom:stability} and reduction~\eqref{eq:axiom:reduction} have first been proved in~\cite{fkmp2013, gantumur2013} with $S(t) = C_{\conststyle stab} \,t$ and $R(t) = C_{\conststyle red} \, t^2$ provided that \eqref{eq:union_of_sons} and~\eqref{eq:reduction:h} are satisfied.
The constants $C_{\conststyle rel}, C_{\conststyle stab} > 0$ depend only on $\Gamma$, $d$, $p$, and uniform $\gamma$-shape regularity of the admissible meshes $\mathcal{T}_H \in \mathbb{T}$, i.e.,
\begin{subequations}\label{eq:Gamma:gamma-shape-regular}
\begin{align}
 \gamma &:= \sup_{\mathcal{T}_H \in \mathbb{T}} \max_{\substack{T,T' \in \mathcal{T}_H\\ T\cap T'\neq\emptyset}} \frac{|T'|}{|T|} < \infty\quad\quad\,\,\, \text{if } d=2,
\intertext{resp.}
 \gamma &:= \sup_{\mathcal{T}_H \in \mathbb{T}} \max_{T \in \mathcal{T}_H} \frac{{\rm diam}(T)}{|T|^{1/(d-1)}} < \infty\quad\text{if } d=3.
\end{align}
\end{subequations}
while $C_{\conststyle red} > 0$ and $0 < \contract{red} < 1$ depend additionally on $\contract{ctr}$.
The {\sl a~priori} convergence follows as in the previous section. Overall, we thus get the subsequent plain convergence result, where we note that for D\"orfler marking~\eqref{eq:marking:doerfler} with sufficiently small marking parameter $0 < \theta \ll 1$, \cite{fkmp2013, gantumur2013} proves even rate-optimal convergence. We also refer to our recent work~\cite{gp2020+}, which proves well-posedness of the residual estimator~\eqref{eq:weaksing:estimator} together with reliability, stability, and reduction for a large class of second-order elliptic PDEs with constant coefficients and the related weakly-singular integral operator $V$ as well as general mesh-refinement strategies.

\begin{proposition}
As long as the mesh-refinement strategy guarantees
regular simplicial triangulations satisfying~\eqref{eq:union_of_sons}, \eqref{eq:reduction:h}, and~\eqref{eq:Gamma:gamma-shape-regular},
and as long as the marking strategy satisfies~\eqref{eq:marking}, Algorithm~\ref{alg:abstract_algorithm} for the weakly-singular integral equation~\eqref{eq:weaksing:weak} driven by the indicators~\eqref{eq:weaksing:estimator} yields convergence
\begin{align}
 C_{\conststyle rel}^{-1} \, \norm{u - u_\ell}{\widetilde{H}^{-1/2}(\Gamma)}
 \le \eta_\ell \to 0
 \quad \text{as } \ell \to \infty.
\end{align}
\end{proposition}

\subsection{Hyper-singular integral equations}

Let $\Omega \subset \mathbb{R}^d$ with $d = 2,3$ be a Lipschitz domain with compact boundary $\partial\Omega$.
On a (relatively) open and connected subset $\Gamma \subseteq \partial\Omega$ and given $f \in H^{-1/2}(\Gamma)$ with $\int_\Gamma f \d{x}=1$
in case of $\Gamma=\partial\Omega$ the hyper-singular integral equation
\begin{align}\label{eq:hypsing:strong}
 (Wu)(x) := \text{p.v.} \int_\Gamma \frac{\partial_x}{\partial\nu(x)}\frac{\partial_y}{\partial\nu(y)}G(x-y) u(y) \d{y} = f(x)
 \quad \text{for all } x \in \Gamma
\end{align}
seeks the unknown integral density $u \in \mathcal{X} := \set{v\in H^{1/2}(\Omega)}{\int_\Gamma v\d{x}=0}$ if $\Gamma=\partial\Omega$ resp.\ $u\in\mathcal{X}:=\widetilde{H}^{1/2}(\Gamma)$ if $\Gamma\subsetneqq\partial\Omega$. Here, $\nu(\cdot)$ denotes the exterior normal vector and $G(\cdot)$ is again the fundamental solution of the Laplacian, i.e., $G(z) = -\frac{1}{2\pi} \log|z|$ for $d = 2$ resp.\ $G(z) = \frac{1}{4\pi} |z|^{-1}$ for $d = 3$.
We note that for $\Gamma = \partial\Omega$, \eqref{eq:hypsing:strong}  is equivalent to the Neumann problem
\begin{align*}
 -\Delta U = 0 \text{ in } \Omega
 \quad \text{subject to}\quad \frac{\partial U}{\partial\nu}  = f \text{ on } \partial\Omega,
\end{align*}
supplemented by the appropriate radiation condition if $\Omega$ is unbounded; see~\cite{mclean}.
The weak formulation reads
\begin{align}\label{eq:hypsing:weak}
 a(u, v) := \int\!\!\!\!\int_{\Gamma\times\Gamma}  G(x-y) \, {\rm curl}_\Gamma u(y) \, {\rm curl}_\Gamma v(x) \d{y} \d{x}
 = \int_\Gamma fv \d{x}
 \quad \text{for all } v \in \mathcal{X},
\end{align}
where ${\rm curl}_\Gamma(\cdot)$ denotes the surface curl (resp.\ the arclength derivative for $d = 2$); see~\cite{mclean}. The
Lax--Milgram lemma yields existence und uniqueness of the solution $u \in \mathcal{X}$.

For a fixed polynomial degree $p \ge 1$ and a regular triangulation $\mathcal{T}_H$ of $\Gamma$ into non-degenerate compact surface simplices, we consider standard boundary element spaces $\mathcal{X}_H = \mathcal{S}^p(\mathcal{T}_H)$ consisting of globally continuous $\mathcal{T}_H$-piecewise polynomials of degree $\le p$ (w.r.t.\ the boundary parametrization). Let $u_H \in \mathcal{X}_H$ be  the corresponding Galerkin solution, i.e.,
\begin{align}
 a(u_H, v_H) = \int_\Gamma f v_H \d{x}
 \quad \text{for all } v_H \in \mathcal{X}_H.
\end{align}
According to the seminal work~\cite{cmps04} and under additional regularity $f \in L^2(\Gamma)$, the local contributions of the residual error estimator read
\begin{align}\label{eq:hypsing:estimator}
 \eta_H(T) := h_T^{1/2} \, \norm{f - Wu_H}{L^2(T)}
 \quad\text{with}\quad
 h_T := |T|^{1/(d-1)},
\end{align}
where $|\cdot|$ is the surface measure. While~\cite{cmps04} proves reliability~\eqref{eq:reliable}, stability~\eqref{eq:axiom:stability} and reduction~\eqref{eq:axiom:reduction} have first been proved in~\cite{gantumur2013,ffkmp15} with $S(t) = C_{\conststyle stab} \,t$ and $R(t) = C_{\conststyle red} \, t^2$ provided that \eqref{eq:union_of_sons} and~\eqref{eq:reduction:h} are satisfied.
The constants $C_{\conststyle rel}, C_{\conststyle stab} > 0$ depend only on $\Gamma$, $d$, $p$, and uniform $\gamma$-shape regularity~\eqref{eq:Gamma:gamma-shape-regular}  of the admissible meshes $\mathcal{T}_H \in \mathbb{T}$,
while $C_{\conststyle red} > 0$ and $0 < \contract{red} < 1$ depend additionally on $\contract{ctr}$.
The {\sl a~priori} convergence follows as in Section~\ref{section:fractional}. Overall, we thus get the subsequent plain convergence result, where we note that for D\"orfler marking~\eqref{eq:marking:doerfler} with sufficiently small marking parameter $0 < \theta \ll 1$, \cite{gantumur2013,ffkmp15} proves even rate-optimal convergence.

\begin{proposition}
As long as the mesh-refinement strategy guarantees regular simplicial triangulations satisfying~\eqref{eq:union_of_sons}, \eqref{eq:reduction:h}, and~\eqref{eq:Gamma:gamma-shape-regular}, and as long as the marking strategy satisfies~\eqref{eq:marking}, Algorithm~\ref{alg:abstract_algorithm} for the hyper-singular integral equation~\eqref{eq:hypsing:weak} driven by the indicators~\eqref{eq:hypsing:estimator} yields convergence
\begin{align}
 C_{\conststyle rel}^{-1} \, \norm{u - u_\ell}{\widetilde{H}^{1/2}(\Gamma)}
 \le \eta_\ell \to 0
 \quad \text{as } \ell \to \infty.
\end{align}
\end{proposition}

\subsection{Nonlinear interface problems}

Let $\Omega \subset \mathbb{R}^d$ with $d = 2,3$ be a bounded Lipschitz domain with compact boundary $\Gamma:=\partial\Omega$ and exterior domain $\Omega^{\rm ext}:=\mathbb{R}^d\setminus\overline\Omega$ such that $\diam(\Omega)<1$ if $d=2$.
Further, let $A:\mathbb{R}^d\to\mathbb{R}^d$ be a Lipschitz continuous and strongly monotone coefficient function in the sense that there exist constants $C_{\rm lip},C_{\rm mon}>0$ such that
\begin{align}
 |Ax-Ay|&\le C_{\rm lip} |x-y| \quad\text{for all }x,y\in\mathbb{R}^d,
 \\
 C_{\rm mon}\norm{\nabla u - \nabla v}{L^2(\Omega)}^2 &\le \int_\Omega (A\nabla u -A \nabla v)\cdot (\nabla u - \nabla v) \d{x}
 \quad\text{for all }u,v\in H^1(\Omega).
\end{align}
For given data $f\in L^2(\Omega)$, $u_D\in H^{1/2}(\Gamma)$, and $\phi_N\in H^{-1/2}(\Gamma)$ with additional compatibility condition
\begin{align}
\int_\Omega f \d{x} + \int_\Gamma \phi_N \d{x} = 0
\quad \text{in case of $d=2$},
\end{align}
we consider the nonlinear interface problem
\begin{subequations}\label{eq:fembem:strong}
\begin{align}
 -\div(A\nabla u) &= f \quad\text{in }\Omega,\\
 -\Delta u^{\rm ext} &= 0 \quad\text{in }\Omega^{\rm ext},\\
 u-u^{\rm ext} &= u_D\quad \text{on }\Gamma,\\
 (A\nabla u - \nabla u^{\rm ext})\cdot \nu &= \phi_N \quad \text{on }\Gamma,\\
 u^{\rm ext} &= \mathcal{O}(|x|^{-1}) \quad \text{as }|x|\to \infty.
\end{align}
\end{subequations}
We seek for a weak solution $(u,u^{\rm ext})\in H^1(\Omega)\times H_{\rm loc}^1(\Omega^{\rm ext})$, where $H_{\rm loc}^1(\Omega^{\rm ext}) = \set{v}{v\in H^1(\omega) \text{ for all open and bounded }\omega\subseteq\Omega^{\rm ext}}$.
There are different ways to equivalently reformulate \eqref{eq:fembem:strong} as FEM-BEM coupling.
To ease presentation, we restrict ourselves to the Bielak--MacCamy coupling~\cite{bm83}, but we stress that Proposition~\ref{prop:fembem:convergence} holds accordingly for  the Johnson--N\'ed\'elec coupling~\cite{jn80} as well as Costabel's symmetric coupling~\cite{costabel88}; see  \cite{affkmp13} for details.
Recalling the single-layer operator $V$ from~\eqref{eq:weaksing:strong} and defining the adjoint double layer operator
\begin{align}
 (K'\phi)(x):=\int_\Gamma \frac{\partial_x}{\partial\nu(x)} G(x,y) \phi(y)  \d{y} \quad\text{for all }\phi\in H^{-1/2}(\Gamma) \text{ and all }  x\in\Gamma,
\end{align}
the variational formulation resulting from the Bielak--MacCamy coupling seeks some ${\bf u} = (u,\phi)\in\mathcal{X}:= H^1(\Omega)\times H^{-1/2}(\Gamma)$ such that
\begin{align}\label{eq:fembem:weak}
 \begin{split}
 a ({\bf u},{\bf v}) &:=
 \int_\Omega (A\nabla u) \cdot \nabla v\d{x} + \int_\Gamma \big((1/2 - K') \phi\big)  v \d{x}  + \int_\Gamma (u - V\phi)\psi \d{x}
 \\
 &=\int_\Omega fv\d{x} + \int_\Gamma\phi_N v \d{x} - \int_\Gamma u_D\psi\d{x} =: F({\bf v})
 \quad \text{for all } {\bf v} = (v,\psi)\in\mathcal{X}.
 \end{split}
\end{align}
According to \cite{affkmp13}, \eqref{eq:fembem:weak} is uniquely solvable provided that $C_{\rm mon}>1/4$.

For a regular triangulation $\mathcal{T}_H$ of $\Omega$ into non-degenerate compact simplices and the induced regular triangulation $\mathcal{E}_H$ of $\Gamma$ into non-degenerate compact surface simplices, we consider globally continuous $\mathcal{T}_H$-piecewise affine functions $\mathcal{S}^1(\mathcal{T}_H)$  to discretize $H^1(\Omega)$, and
 $\mathcal{E}_H$-piecewise constant functions $\mathcal{P}^0(\mathcal{E}_H)$ to discretize $H^{-1/2}(\Gamma)$,
 i.e., $\mathcal{X}_H=\mathcal{S}^1(\mathcal{T}_H)\times\mathcal P^0(\mathcal{E}_H)$.
  Let ${\bf u}_H \in \mathcal{X}_H$ be  the corresponding Galerkin solution, i.e.,
\begin{align}
 a({\bf u}_H, {\bf v}_H) =F({\bf v}_H)
 \quad \text{for all } {\bf v}_H \in \mathcal{X}_H.
\end{align}
According to~\cite{affkmp13}, the local contributions of the residual error estimator read
\begin{align}\label{eq:fembem:estimator}
 \begin{split}
 &\eta_H(T)^2 :=
 h_T^2 \norm{f}{L^2(T)}^2 + h_T \Big(\norm{[(A\nabla u_H)\cdot \nu]}{L^2(\partial T\cap \Omega)}^2
 \\
 &+ \norm{\phi_N+(K'-1/2)\phi_H-(A\nabla u_H)\cdot\nu}{L^2(\partial T\cap\Gamma)}^2 + \norm{\nabla_\Gamma (u_H-u_D-V\phi_H)}{L^2(\partial T\cap\Gamma)}^2\Big)
 \end{split}
\end{align}
with the surface gradient $\nabla_\Gamma(\cdot)$ and the mesh-size $h_T:=|T|^{1/d}$ for all $T\in\mathcal{T}_H$.
Indeed, \cite{affkmp13} proves reliability~\eqref{eq:reliable}, stability~\eqref{eq:axiom:stability}, and reduction~\eqref{eq:axiom:reduction} (where the terms $u_H,u_h$ are replaced by ${\bf u}_H,{\bf u}_h$)  with $S(t) = C_{\conststyle stab} \,t$ and $R(t) = C_{\conststyle red} \, t^2$ provided that \eqref{eq:union_of_sons} and~\eqref{eq:reduction:h} are satisfied.
The constants $C_{\conststyle rel}, C_{\conststyle stab} > 0$ depend only on $\Gamma$, $d$, $C_{\rm lip}$, and uniform $\gamma$-shape regularity~\eqref{eq:gamma-shape-regular}  of the admissible meshes $\mathcal{T}_H \in \mathbb{T}$,
while $C_{\conststyle red} > 0$ and $0 < \contract{red} < 1$ depend additionally on $\contract{ctr}$.
The {\sl a~priori} convergence follows as in Section~\ref{section:fractional}, where the required C\'ea lemma is given in~\cite[Corollary~12]{affkmp13}. Overall, we thus get the subsequent plain convergence result.

\begin{proposition}\label{prop:fembem:convergence}
As long as the mesh-refinement strategy guarantees regular simplicial triangulations satisfying~\eqref{eq:union_of_sons}, \eqref{eq:reduction:h}, and~\eqref{eq:gamma-shape-regular}, and as long as the marking strategy satisfies~\eqref{eq:marking}, Algorithm~\ref{alg:abstract_algorithm} for the Bielak--MacCamy coupling~\eqref{eq:fembem:weak} driven by the indicators~\eqref{eq:fembem:estimator} yields convergence
\begin{align}
 C_{\conststyle rel}^{-1} \, \norm{{\bf u} - {\bf u}_\ell}{H^1(\Omega)\times H^{-1/2}(\Gamma)}
 \le \eta_\ell \to 0
 \quad \text{as } \ell \to \infty.
\end{align}
\end{proposition}

\section{\@startsection{section}{1}%
\z@{.7\linespacing\@plus\linespacing}{.5\linespacing}%
{\normalsize\scshape\bfseries\centering}}
\makeatother

\makeatletter
\renewcommand{\@secnumfont}{\bfseries}
\makeatother

\newcounter{statement}
\newenvironment{statement}[2][!]{%
\vskip2mm
\noindent%
\refstepcounter{statement}%
\bf#2~\thestatement%
\ifthenelse{\equal{#1}{!}}{.\ }{~(#1).\ }%
\it%
}{%
\vskip1mm
}

\newenvironment{theorem}[1][!]{\begin{statement}[#1]{Theorem}}{\end{statement}}
\newenvironment{lemma}[1][!]{\begin{statement}[#1]{Lemma}}{\end{statement}}
\newenvironment{definition}[1][!]{\begin{statement}[#1]{Definition}}{\end{statement}}
\newenvironment{proposition}[1][!]{\begin{statement}[#1]{Proposition}}{\end{statement}}
\newenvironment{corollary}[1][!]{\begin{statement}[#1]{Corollary}}{\end{statement}}
\newenvironment{remark}[1][!]{\begin{statement}[#1]{Remark}}{\end{statement}}
\newenvironment{algorithm}[1][!]{\begin{statement}[#1]{Algorithm}}{\end{statement}}

\defH{H}
\defh{h}
\def\mathcal{X}{\mathcal{X}}
\def\mathcal{U}{\mathcal{U}}
\def\mathbb{T}{\mathbb{T}}
\def\const#1{C_{\conststyle#1}}
\def\contract#1{q_{\conststyle#1}}
\def\reff#1#2{\stackrel{\eqref{#1}}{#2}}
\def\phantomreff#1#2{\stackrel{\phantom{\eqref{#1}}}{#2}}

\begin{document}

\begin{abstract}
We consider $h$-adaptive algorithms in the context of the finite element method (FEM) and the boundary element method (BEM). Under quite general assumptions on the building blocks SOLVE, ESTIMATE, MARK, and REFINE of such algorithms, we prove plain convergence in the sense that the adaptive algorithm drives the underlying {\sl a~posteriori} error estimator to zero. Unlike available results in the literature, our analysis avoids the use of any reliability and efficiency estimate, but only relies on structural properties of the estimator, namely stability on non-refined elements and reduction on refined elements. In particular, the new framework thus covers also problems involving non-local operators like the fractional Laplacian or boundary integral equations, where (discrete) efficiency is (currently) not available.
\end{abstract}

\maketitle
\thispagestyle{fancy}

\def\mathcal{S}{\mathcal{S}}
\def{\color{red}XXX}{{\color{red}XXX}}
\section{Introduction}

\noindent
{\sl A~posteriori} error estimation and related adaptive mesh-refinement via the loop
\begin{align}\label{eq:semr}
\bf SOLVE
\quad \longrightarrow \quad
ESTIMATE
\quad \longrightarrow \quad
MARK
\quad \longrightarrow \quad
REFINE
\end{align}
are standard tools in modern scientific computing. Over the last decade, the mathematical understanding has matured.
Convergence with optimal algebraic rates is mathematically guaranteed for a reasonable class of elliptic model problems and standard discretizations; we refer to the  works~\cite{doerfler1996,mns2000,bdd2004,stevenson2007,ckns2008,cn2012,ffp2014} for some important steps as well as to the state-of-the-art review~\cite{axioms}. However, all these works employ the so-called D\"orfler marking strategy proposed in~\cite{doerfler1996} to single out elements for refinement. Moreover, for the 2D Poisson problem, it has recently been shown that a modified maximum criterion does not only lead to optimal convergence rates, but even leads to instance optimal meshes~\cite{dks2016,ks2016,ip2020+}. As the focus comes to other marking strategies, only plain convergence results are known and the essential works are~\cite{msv2008,siebert2011}.

To outline the results of~\cite{msv2008,siebert2011} and the contributions of the present work, let us fix some notation.
Let $\mathcal{X}$ be a normed space which is linked to some domain (or manifold) $\Omega \subset \mathbb{R}^d$, $d \ge 1$. Let $u \in \mathcal{X}$  be the sought (unknown) solution. Suppose that the discrete subspaces $\mathcal{X}_\ell \subset \mathcal{X}$ are linked to some mesh $\mathcal{T}_\ell$ of $\Omega$ consisting of compact subdomains of $\Omega$. Let $u_\ell \in \mathcal{X}_\ell$ be a \emph{computable} discrete approximation of $u$. Finally, let $\eta_\ell^2 = \sum_{T \in \mathcal{T}_\ell} \eta_\ell(T)^2$ be a \emph{computable} error estimator such that $\eta_\ell(T)$ measures, at least heuristically, the error $u - u_\ell$ on $T \in \mathcal{T}_\ell$. We suppose that the sequence of meshes $(\mathcal{T}_\ell)_{\ell \in \mathbb{N}_0}$ is generated by the adaptive loop~\eqref{eq:semr}. In such a setting, it has already been observed in the seminal work~\cite{bv1984} that nestedness $\mathcal{X}_\ell \subseteq \mathcal{X}_{\ell+1}$ of the discrete spaces together with a C\'ea-type quasi-optimality proves the so-called {\sl a~priori} convergence of adaptive schemes, i.e., there always exists a limit $u_\infty \in \mathcal{X}$ such that
\begin{align}\label{assumption:limit_u}
 \norm{u_\infty - u_\ell}{\mathcal{X}} \to 0
 \quad \text{as } \ell \to \infty.
\end{align}
However, it remains to prove that also $u = u_\infty$.

To explain the abstract notation, let us consider the 2D Poisson model problem: In this case,
 $\Omega \subset \mathbb{R}^2$ is a polygonal Lipschitz domain, $f \in L^2(\Omega)$ is some given load, $u \in \mathcal{X} = H^1_0(\Omega)$ solves the 2D Poisson model problem $-\Delta u = f$ in $\Omega$ subject to homogeneous boundary conditions $u = 0$ on $\partial\Omega$, the meshes $\mathcal{T}_\ell$ are conforming triangulations of $\Omega$ into compact triangles $T \in \mathcal{T}_\ell$, and $u_\ell \in \mathcal{X}_\ell = \set{v_\ell \in H^1_0(\Omega)}{v_\ell|_T \text{ is affine for all } T \in \mathcal{T}_\ell}$ is the conforming first-order finite element approximation of $u$, which solves
\begin{align}\label{eq:poisson:fem}
 \int_\Omega \nabla u_\ell \cdot \nabla v_\ell \d{x} = \int_\Omega f v_\ell \d{x}
 \quad \text{for all } v_\ell \in \mathcal{X}_\ell
\end{align}
For this problem, the classical residual error estimator reads
\begin{align}\label{eq:poisson:estimator}
 \eta_\ell^2 = \sum_{T \in \mathcal{T}_\ell} \eta_\ell(T)^2
 \quad \text{with} \quad
 \eta_\ell(T)^2 = h_T^2 \, \norm{f}{L^2(T)}^2 + h_T \, \norm{\llbracket\partial_n u_\ell\rrbracket}{L^2(\partial T \cap \Omega)}^2,
\end{align}
where $\llbracket\cdot\rrbracket$ denotes the jump across interior edges and $h_T = |T|^{1/2}$ denotes the local mesh-size; see, e.g., the monographs~\cite{ao00,verfuerth}.

While~\cite{msv2008} formally focusses on conforming Petrov--Galerkin discretizations in the setting of Ladyshenskaja--Babuska--Brezzi (LBB), the actual analysis is more general: Besides some assumptions on the locality of the norms of the involved function spaces, there are no assumptions on how $u$ or $u_\ell$ are computed. The crucial assumptions in~\cite{msv2008} are \emph{local efficiency}
\begin{align}\label{eq:local_efficiency}
 C_{\rm eff}^{-1} \, \eta_\ell(T)
 \le \norm{u - u_\ell}{\mathcal{X}(\Omega_\ell(T))} + {\rm osc}_\ell(\Omega_\ell(T))
 \quad \text{for all } T \in \mathcal{T}_\ell
\end{align}
as well as \emph{discrete local efficiency} on marked elements
\begin{align}\label{eq:discrete_local_efficiency}
 C_{\rm eff}^{-1} \, \eta_\ell(T)
 \le \norm{u_{\ell+1} - u_\ell}{\mathcal{X}(\Omega_\ell(T))} + {\rm osc}_\ell(\Omega_\ell(T))
 \quad \text{for all } T \in \mathcal{M}_\ell,
\end{align}
where $\Omega_\ell(T) = \bigcup\set{T' \in \mathcal{T}_\ell}{T' \cap T \neq \emptyset}$ is the patch of $T$ and ${\rm osc}_\ell$ are some data oscillation terms. It is known that the latter assumption requires (at least) stronger local refinement, e.g., the local bisec5 refinement of marked elements in 2D to ensure the interior node property; see, e.g.,~\cite{mns2000}. The main result of~\cite{msv2008} proves that under these assumptions and for quite general marking strategies (see~\eqref{eq:marking} below), the adaptive algorithm ensures that {\sl a~priori} convergence~\eqref{assumption:limit_u} already implies \emph{estimator convergence}
\begin{align}\label{eq:estimator_convergence}
 \eta_\ell \to 0
 \quad \text{as } \ell \to \infty.
\end{align}
Provided that the error estimator $\eta_\ell$ additionally satisfies reliability, i.e.,
\begin{align}\label{eq:reliability}
 \norm{u - u_\ell}{\mathcal{X}} \le \const{rel} \, \eta_\ell,
\end{align}
this proves that $\norm{u - u_\ell}{\mathcal{X}} \to 0$ as $\ell \to \infty$. In explicit terms, the main result of~\cite{msv2008} reads as follows: If the discrete solutions $u_\ell \in \mathcal{X}_\ell$ converge~\eqref{assumption:limit_u}, then they converge indeed to the correct limit $u = u_\infty$ --- provided that the error estimator satisfies~\eqref{eq:local_efficiency}--\eqref{eq:discrete_local_efficiency} and~\eqref{eq:reliability}.

Conceptually, it is remarkable that the convergence proof of~\cite{msv2008} exploits lower error bounds, although the mesh-refinement is driven by the error estimator only. The work~\cite{siebert2011} thus aimed to prove convergence without using (discrete) lower bounds. This, however, comes at the cost that, first, the analysis exploits the problem setting (and is restricted to Petrov--Galerkin discretizations of operator equations $Bu = F$) and, second, the ana\-ly\-sis relies on some strengthened reliability estimate (formulated in terms of the residual), which implies~\eqref{eq:reliability}. The main result of~\cite{siebert2011} then states that under these assumptions and for quite general marking strategies (see~\eqref{eq:marking2} below), the adaptive algorithm ensures that {\sl a~priori} convergence~\eqref{assumption:limit_u} already implies \emph{error convergence}
\begin{align}\label{eq:error_convergence}
 \norm{u - u_\ell}{\mathcal{X}} \to 0
 \quad \text{as } \ell \to \infty.
\end{align}
In particular, the new proof of~\cite{siebert2011} avoids the discrete local efficiency~\eqref{eq:discrete_local_efficiency}.
Surprisingly, however, estimator convergence~\eqref{eq:estimator_convergence} cannot be proved under the assumptions of~\cite{siebert2011} but requires that the error estimator $\eta_\ell$ is also locally efficient~\eqref{eq:local_efficiency}.

One advantage of the results of~\cite{msv2008,siebert2011} is that they apply to  many different {\sl a~posteriori} error estimators. In particular, it has recently been shown in~\cite{fp2020+,gs2020+} that the assumptions of~\cite{siebert2011} are, in particular, satisfied for a wide range of model problems discretized by least squares finite element methods, where adaptivity is driven by the built-in least-squares functional, including even a least-squares space-time discretization of the heat equation. On the other hand,~\cite{siebert2011} excludes adaptive schemes for variational inequalities, and both works~\cite{msv2008,siebert2011} need local efficiency of the error estimator, which does not appear to be available for non-local operators, e.g., finite element methods for the fractional Laplacian (see, e.g.,~\cite{fmp2019+}) or boundary element methods for elliptic integral equations (see, e.g.,~\cite{fkmp2013,gantumur2013}).

With the latter observations, the current paper comes into play. We provide a new proof for plain convergence of adaptive algorithms, which does \emph{neither} involve reliability~\eqref{eq:reliability} \emph{nor} any kind of (global or local) efficiency~\eqref{eq:local_efficiency}--\eqref{eq:discrete_local_efficiency}. Instead, we exploit that the local contributions of many residual error estimators are weighted by the local mesh-size (cf.~\eqref{eq:poisson:estimator} for the Poisson model problem). With scaling arguments, one usually obtains \emph{reduction on refined elements}
\begin{align}\label{eq:intro:reduction}
 \eta_{\ell+n}(\mathcal{T}_{\ell+n} \backslash \mathcal{T}_\ell)^2
 \le q \, \eta_{\ell}(\mathcal{T}_\ell \backslash \mathcal{T}_{\ell+n})^2
 + C \, \norm{u_{\ell+n} - u_\ell}{\mathcal{X}}^2
 \quad \text{for all } \ell, n \in \mathbb{N}_0,
\end{align}
with generic constants $0 < q < 1$ and $C > 0$. Note that $\mathcal{T}_\ell \backslash \mathcal{T}_{\ell+n}$ corresponds to the elements that are going to be refined, while $\mathcal{T}_{\ell+n} \backslash \mathcal{T}_\ell$ corresponds to the generated children. Moreover, there usually holds \emph{stability on non-refined elements}
\begin{align}\label{eq:intro:stability}
 | \eta_{\ell+n}(\mathcal{T}_{\ell+n} \cap \mathcal{T}_\ell) - \eta_\ell(\mathcal{T}_{\ell+n} \cap \mathcal{T}_\ell) |
 \le C \, \norm{u_{\ell+n} - u_\ell}{\mathcal{X}}
 \quad \text{for all } \ell,n \in \mathbb{N}_0.
\end{align}
We stress that~\eqref{eq:intro:reduction}--\eqref{eq:intro:stability} play also a fundamental role in the contemporary proofs of optimal convergence rates for adaptive algorithms; see~\cite{axioms}. The main result of the present work (Theorem~\ref{theoerem:estconv:new}) shows that, together with the same marking criterion as in~\cite{msv2008}, the structural properties~\eqref{eq:intro:reduction}--\eqref{eq:intro:stability} suffice to show that {\sl a~priori} convergence~\eqref{assumption:limit_u} yields estimator convergence~\eqref{eq:estimator_convergence}. Clearly, reliability~\eqref{eq:reliability} is then finally required to conclude error convergence~\eqref{eq:error_convergence}.

\bigskip
\textbf{Outline.} The remainder of this work is organized as follows: In Section~\ref{section:abstract}, we provide a formal statement of the adaptive algorithm (Algorithm~\ref{alg:abstract_algorithm}) as well as precise assumptions on its four modules from~\eqref{eq:semr}. The new plain convergence result (Theorem~\ref{theoerem:estconv:new}, Theorem~\ref{theoerem2:estconv:new}) is stated and proved in Section~\ref{section:main}, before we give some examples which do not fit the framework of~\cite{msv2008,siebert2011}, but are covered by the current analysis.

\section{Abstract adaptive algorithm}
\label{section:abstract}

\subsection{Mesh refinement}\label{section:mesh-refinement}

Let $\Omega \subset \mathbb{R}^d$ be a bounded domain (or a manifold in $\mathbb{R}^d$) with positive measure $|\Omega| > 0$. We say that $\mathcal{T}_H$ is a mesh (of $\Omega)$, if
\begin{itemize}
\item $\mathcal{T}_H$ is a finite set of compact sets $T \in \mathcal{T}_H$ with positive measure $|T| > 0$;
\item for all $T, T' \in \mathcal{T}_H$ with $T \neq T'$, it holds that $|T \cap T'| = 0$;
\item $\mathcal{T}_H$ is a covering of $\overline\Omega$, i.e., $\overline\Omega = \bigcup_{T \in \mathcal{T}_H} T$.
\end{itemize}
Let ${\rm refine}(\cdot)$ be a fixed refinement strategy, i.e., for each mesh $\mathcal{T}_H$ and a set of marked elements $\mathcal{M}_H \subseteq \mathcal{T}_H$, the refinement strategy returns a refined mesh $\mathcal{T}_h := {\rm refine}(\mathcal{T}_H,\mathcal{M}_H)$ such that, first, at least the marked elements are refined (i.e., $\mathcal{M}_H \subseteq \mathcal{T}_H \backslash \mathcal{T}_h$) and, second,
parents $T \in \mathcal{T}_H$ are the union of their children, i.e.,
\begin{align}\label{eq:union_of_sons}
 T = \bigcup \set{T' \in \mathcal{T}_h}{T' \subseteq T}
 \quad \text{for all } T \in \mathcal{T}_H.
\end{align}
For a mesh $\mathcal{T}_H$, let $\mathbb{T}(\mathcal{T}_H)$ denote the set of all possible refinements of $\mathcal{T}_H$ (as determined by the refinement strategy ${\rm refine}(\cdot)$), i.e., for any $\mathcal{T}_h \in \mathbb{T}(\mathcal{T}_H)$, there exists $n \in \mathbb{N}_0$ and $\mathcal{T}'_0, \dots, \mathcal{T}'_n$ such that $\mathcal{T}'_0 = \mathcal{T}_H$, $\mathcal{T}'_{j+1} = {\rm refine}(\mathcal{T}'_j,\mathcal{M}'_j)$ for all $j = 0, \dots, n-1$ and appropriate $\mathcal{M}'_j \subseteq \mathcal{T}'_j$, and $\mathcal{T}'_n = \mathcal{T}_h$.
Finally, we suppose that we are given a fixed initial mesh $\mathcal{T}_0$ so that it makes sense to call $\mathbb{T} := \mathbb{T}(\mathcal{T}_0)$ the \emph{set of all admissible meshes.}

\subsection{Continuous and discrete setting}

Let $\mathcal{X}$ be a normed space (related to $\Omega$) and $u \in \mathcal{X}$ be the (unknown) exact solution. For each mesh $\mathcal{T}_H$, let $\mathcal{X}_H \subseteq \mathcal{X}$ be an associated discrete subspace and $u_H \in \mathcal{X}_H$ be the corresponding (computable) discrete solution.

\subsection{Error estimator}

For each mesh $\mathcal{T}_H$ and all $T \in \mathcal{T}_H$, let $\eta_H(T) \ge 0$ be a \emph{computable} quantity which is usually called \emph{refinement indicator}. At least heuristically, $\eta_H(T)$ measures the error $u - u_H$ on the element $T$. We abbreviate
\begin{align}
 \eta_H := \eta_H(\mathcal{T}_H),
 \quad \text{where} \quad
 \eta_H(\mathcal{U}_H) := \Big( \sum_{T \in \mathcal{T}_H} \eta_H(T)^2 \Big)^{1/2}
 \quad \text{for all } \mathcal{U}_H \subseteq \mathcal{T}_H.
\end{align}
We note that $\eta_H$ is usually referred to as \emph{error estimator}.

\subsection{Adaptive algorithm}

Starting from the given initial mesh $\mathcal{T}_0$, we consider the standard adaptive loop~\eqref{eq:semr} in the following algorithmic form:

\begin{algorithm}
\label{alg:abstract_algorithm}
For each $\ell=0,1,2,\dots$, iterate the following steps~{\rm(i)--(iv)}:
\begin{itemize
\item[\rm(i)] {\bf SOLVE:} Compute the discrete solution $u_\ell \in \mathcal{X}_\ell$.
\item[\rm(ii)] {\bf ESTIMATE:} Compute  refinement indicators $\eta_\ell(T)$ for all elements ${T}\in\mathcal{T}_\ell$.
\item[\rm(iii)] {\bf MARK:} Determine a set of marked elements $\mathcal{M}_\ell\subseteq\mathcal{T}_\ell$.
\item[\rm(iv)] {\bf REFINE:} Generate the refined mesh $\mathcal{T}_{\ell+1}:={\rm refine}(\mathcal{T}_\ell,\mathcal{M}_\ell)$.
\end{itemize}
\textbf{Output:} Refined meshes $\mathcal{T}_\ell$, corresponding exact discrete solutions $u_\ell$, and
error estimators $\eta_\ell$ for all $\ell \in \mathbb{N}_0$.\qed
\end{algorithm}

To analyze Algorithm~\ref{alg:abstract_algorithm}, it remains to specify further assumptions on its four modules: As far as {\bf SOLVE} is concerned, we shall only assume {\sl a~priori} convergence~\eqref{assumption:limit_u}.
While this assumption is guaranteed for many problems (see, e.g.,~\cite{msv2008,siebert2011} for problems in the framework of the LBB theory as well as the seminal work~\cite{bv1984} for problems in the Lax--Milgram setting), we stress that, at this point, it is still mathematically unclear whether there holds $u_\infty = u$ or not.

As far as {\bf MARK} is concerned, let $M \colon \mathbb{R}_{\ge0} \to \mathbb{R}_{\ge0}$ be continuous at $0$ with $M(0)=0$ and suppose that the sets $\mathcal{M}_\ell \subseteq \mathcal{T}_\ell$ satisfy the following property from~\cite{msv2008}:
\begin{align}\label{eq:marking}
 \max_{T\in\mathcal{T}_\ell\setminus\mathcal{M}_\ell} \eta_\ell(T) \le M\big(\eta_\ell(\mathcal{M}_\ell)\big).
\end{align}
We note that the latter assumption is weaker than the following assumption from~\cite{siebert2011}:
\begin{align}\label{eq:marking2}
 \max_{T\in\mathcal{T}_\ell\setminus\mathcal{M}_\ell} \eta_\ell(T) \le M\big(\max_{T \in \mathcal{M}_\ell}\eta_\ell(T)\big).
\end{align}
Clearly, the marking criteria~\eqref{eq:marking}--\eqref{eq:marking2} are satisfied with $M(t) = t$ as soon as $\mathcal{M}_\ell$ contains one element with maximal indicator, i.e., there exists $T \in \mathcal{M}_\ell$ such that $\eta_\ell(T) = \max_{T' \in \mathcal{T}_\ell} \eta_\ell(T')$. For instance, this is the case for
\begin{itemize}
\item the \emph{maximum criterion} for some fixed $0 \le \theta \le 1$, where
\begin{align}\label{eq:marking:maximum}
 \mathcal{M}_\ell:=\set{T\in\mathcal{T}_\ell}{\eta_\ell(T) \ge (1-\theta) \max_{T'\in\mathcal{T}_\ell} \eta_\ell(T')}
\end{align}
\item the \emph{equidistribution criterion} for fixed $0 \le \theta \le 1$, where
\begin{align}\label{eq:marking:equilibration}
 \mathcal{M}_\ell := \set{T\in\mathcal{T}_\ell}{\eta_\ell(T)\ge (1-\theta)\, \eta_\ell/\#\mathcal{T}_\ell}
\end{align}
\end{itemize}
Finally, let us consider the \emph{D\"orfler criterion} for some fixed $0<\theta\le1$ , i.e.,
\begin{align}\label{eq:marking:doerfler}
 \theta\,\eta_\ell^2 \le \eta_\ell(\mathcal{M}_\ell)^2.
\end{align}
While~\eqref{eq:marking2} cannot be satisfied in general,~\eqref{eq:marking} holds with $M(t):=\sqrt{(1-\theta)\,\theta^{-1}t}$.
To see this, let $T\in\mathcal{T}_\ell\setminus\mathcal{M}_\ell$ and note that
\begin{align*}
 \eta_\ell(T)^2
 &\le \eta_\ell(\mathcal{T}_\ell\setminus\mathcal{M}_\ell)^2

 = \eta_\ell^2 - \eta_\ell(\mathcal{M}_\ell)^2 \le (1-\theta)\eta_\ell^2 \le (1-\theta)\,\theta^{-1} \eta_\ell(\mathcal{M}_\ell)^2.
\end{align*}
However, if the set $\mathcal{M}_\ell$ is constructed via sorting of the indicators, then
\begin{align}\label{eq:marking:doerfler:max}
 \max_{T'\in\mathcal{T}_\ell\setminus\mathcal{M}_\ell} \eta_\ell(T')
 \le \min_{T'\in\mathcal{M}_\ell} \eta_\ell(T');
\end{align}
see \cite{pp2020} for different algorithms which generate $\mathcal{M}_\ell \subseteq \mathcal{T}_\ell$ satisfying the D\"orfler criterion~\eqref{eq:marking:doerfler} together with~\eqref{eq:marking:doerfler:max}. In the latter case,~\eqref{eq:marking}--\eqref{eq:marking2} hold again with $M(t):=t$.

\section{A new plain convergence result}
\label{section:main}

Unlike~\cite{msv2008,siebert2011}, we only require the following two structural properties of the error estimator for all $\mathcal{T}_H \in \mathbb{T}$ and all refinements $\mathcal{T}_h \in \mathbb{T}(\mathcal{T}_H)$, where $S, R \colon \mathbb{R}_{\ge0} \to \mathbb{R}_{\ge0}$ are functions which are continuous at $0$ with $R(0) = 0 = S(0)$ and $0 < \contract{red} < 1$:
\begin{itemize}
\item {\bf stability on non-refined elements}, i.e.,
\begin{align}\label{eq:axiom:stability}
 \eta_h(\mathcal{T}_h \cap \mathcal{T}_H)
 \le \eta_H(\mathcal{T}_h \cap \mathcal{T}_H)
 + S(\norm{u_h - u_H}{\mathcal{X}});
\end{align}
\item {\bf reduction on refined elements}, i.e.,
\begin{align}\label{eq:axiom:reduction}
 \eta_h(\mathcal{T}_h \backslash \mathcal{T}_H)^2
 \le
 \contract{red} \, \eta_H(\mathcal{T}_H \backslash \mathcal{T}_h)^2
 + R(\norm{u_h - u_H}{\mathcal{X}}).
\end{align}
\end{itemize}
We note that, \eqref{eq:axiom:stability}--\eqref{eq:axiom:reduction} are implicitly first found in the proof of~\cite[Corollary~3.4]{ckns2008}, but already seem to go back to~\cite{dk2008} (used there for the oscillations).
In practice, the reduction~\eqref{eq:axiom:reduction} can only be proved if the local contributions $\eta_H(T)$ of the error estimator are weighted by (some positive power of) the local mesh-size $h_T$. In the later examples in Section~\ref{section:examples}, it holds that $S(t) \sim t$ and $R(t) \sim t^2$.

Under the structural assumptions~\eqref{eq:axiom:stability}--\eqref{eq:axiom:reduction} on the estimator, the following theorem already proves that Algorithm~\ref{alg:abstract_algorithm} leads to estimator convergence. We stress that neither the reliability estimate~\eqref{eq:reliability} nor any (global or even local) efficiency estimate (e.g.,~\eqref{eq:local_efficiency}--\eqref{eq:discrete_local_efficiency}) is required.

\begin{theorem}\label{theoerem:estconv:new}
Suppose the properties~\eqref{eq:axiom:stability}--\eqref{eq:axiom:reduction} of the estimator and that refinement ensures that each parent is the union of its children~\eqref{eq:union_of_sons}. Consider the output of Algorithm~\ref{alg:abstract_algorithm}
with the marking strategy~\eqref{eq:marking}.
Then, {\sl a~priori} convergence~\eqref{assumption:limit_u} implies estimator convergence
\begin{align}\label{eq:theoerem:estconv:new}
 \eta_\ell \to 0 \quad \text{as } \ell \to \infty.
\end{align}
\end{theorem}

The proof of Theorem~\ref{theoerem:estconv:new} employs the following elementary result, whose simple proof is included for the convenience of the reader.

\begin{lemma}\label{axioms:lemma:cor:estimator_convergence}
Let $(a_\ell)_{\ell \in \mathbb{N}_0}$ be a sequence with $a_\ell \ge 0$ for all $\ell \in \mathbb{N}_0$.
Suppose that there exists $0 < \rho < 1$ and a sequence $(b_\ell)_{\ell \in \mathbb{N}_0}$ with $b_\ell \to 0$ as $\ell \to \infty$ such that
\begin{align}
 a_{\ell+1} \le \rho a_\ell + b_\ell
 \quad \text{for all } \ell \in \mathbb{N}_0.
\end{align}
Then, if follows that $a_\ell \to 0$ as $\ell \to \infty$.
\end{lemma}

\begin{proof}
With the convergence of $(b_\ell)_{\ell \in \mathbb{N}_0}$, we note that
\begin{align*}
0 \le \limsup_{\ell\to\infty} a_\ell
= \limsup_{\ell\to\infty} a_{\ell+1}
&\le \limsup_{\ell\to\infty} \, ( \rho a_\ell + b_\ell ) = \rho \, \limsup_{\ell\to\infty}a_\ell.
\end{align*}
Thus, it only remains to show that $\limsup_{\ell\to\infty}a_\ell < \infty$ to conclude that $0=\liminf_{\ell\to\infty} a_\ell = \limsup_{\ell\to\infty}a_\ell$ and hence $\lim_{\ell\to\infty} a_\ell = 0$.
Indeed, induction on $\ell$ proves  that
\begin{align*}
0 \le a_{\ell}
&\le \rho^{\ell} a_0 + \sum_{j=0}^{\ell-1} \rho^{\ell-1-j} b_j
\quad \text{for all $\ell\in\mathbb{N}_0$.}
\end{align*}
Since $(b_\ell)_{\ell \in \mathbb{N}_0}$ is uniformly bounded, the geometric series yields that $\sup_{\ell \in \mathbb{N}} a_\ell < \infty$.
In particular, we thus see that $\limsup_{\ell\to\infty}a_\ell < \infty$ and conclude the proof.
\end{proof}

\def\varepsilon{\varepsilon}
\begin{proof}[Proof of Theorem~\ref{theoerem:estconv:new}]
The proof is split into five steps.

{\bf Step 1:} We prove that $(\mathcal{T}_\ell)_{\ell \in \mathbb{N}_0}$ admits a subsequence $(\mathcal{T}_{\ell_k})_{k \in \mathbb{N}_0}$ such that
$\eta_{\ell_k}(\mathcal{T}_{\ell_k} \backslash \mathcal{T}_{\ell_{k+1}}) \to 0$ as $k \to \infty$.
Let $\mathcal{T}_\infty = \bigcup_{\ell' \in \mathbb{N}_0} \bigcap_{\ell \ge \ell'} \mathcal{T}_\ell$ be the set of all elements, which remain unrefined after some (arbitrary) step $\ell'$. We exploit~\eqref{eq:union_of_sons} and choose a subsequence $(\mathcal{T}_{\ell_k})_{k \in \mathbb{N}_0}$ of $(\mathcal{T}_\ell)_{\ell \in \mathbb{N}_0}$ such that
\begin{align}\label{eq10:theoerem:estconv:new}
 \mathcal{T}_{\ell_{k+1}} \cap \mathcal{T}_{\ell_k} = \mathcal{T}_{\ell_k} \cap \mathcal{T}_\infty
 \quad \text{for all } k \in \mathbb{N}_0,
\end{align}
i.e., only elements $T \in \mathcal{T}_{\ell_k} \cap \mathcal{T}_\infty$ remain unrefined if we pass from $\mathcal{T}_{\ell_k}$ to $\mathcal{T}_{\ell_{k+1}}$.
Note that the choice of $(\mathcal{T}_{\ell_k})_{k \in \mathbb{N}_0}$ guarantees the inclusion
\begin{align*}
 \mathcal{T}_{\ell_{k+1}} \cap \mathcal{T}_{\ell_k}
 \reff{eq10:theoerem:estconv:new}= \mathcal{T}_{\ell_k} \cap \mathcal{T}_\infty
 \subseteq \mathcal{T}_{\ell_{k+1}} \cap \mathcal{T}_\infty
 \reff{eq10:theoerem:estconv:new}= \mathcal{T}_{\ell_{k+2}} \cap \mathcal{T}_{\ell_{k+1}}
 \end{align*}
 and hence
\begin{align*}
 \mathcal{T}_{\ell_{k+1}} \backslash \mathcal{T}_{\ell_{k+2}}
 = \mathcal{T}_{\ell_{k+1}} \backslash [ \mathcal{T}_{\ell_{k+2}} \cap \mathcal{T}_{\ell_{k+1}} ]
 \subseteq \mathcal{T}_{\ell_{k+1}} \backslash [ \mathcal{T}_{\ell_{k+1}} \cap \mathcal{T}_{\ell_k} ]
 = \mathcal{T}_{\ell_{k+1}} \backslash \mathcal{T}_{\ell_k}.
\end{align*}
With this and reduction~\eqref{eq:axiom:reduction}, we infer that
\begin{align*}
 \eta_{\ell_{k+1}}(\mathcal{T}_{\ell_{k+1}} \backslash \mathcal{T}_{\ell_{k+2}})^2
 &\le \eta_{\ell_{k+1}}(\mathcal{T}_{\ell_{k+1}} \backslash \mathcal{T}_{\ell_k})^2

 \reff{eq:axiom:reduction}\le \contract{red} \, \eta_{\ell_k}(\mathcal{T}_{\ell_k} \backslash \mathcal{T}_{\ell_{k+1}})^2
 + R(\norm{u_{\ell_{k+1}} - u_{\ell_k}}{\mathcal{X}}).
\end{align*}
With $a_k = \eta_{\ell_k}(\mathcal{T}_{\ell_k} \backslash \mathcal{T}_{\ell_{k+1}})^2$, $\rho := \contract{red}$, and $b_k := R(\norm{u_{\ell_{k+1}} - u_{\ell_k}}{\mathcal{X}})$, {\sl a~priori} convergence~\eqref{assumption:limit_u} proves that
\begin{align*}
 0 \le a_{k+1} \le \rho \, a_k + b_k
 \quad \text{for all } k \in \mathbb{N}_0
 \quad \text{with } \lim_{k \to \infty} b_k = 0.
\end{align*}
By use of Lemma~\ref{axioms:lemma:cor:estimator_convergence}, we conclude that $a_k = \eta_{\ell_k}(\mathcal{T}_{\ell_k} \backslash \mathcal{T}_{\ell_{k+1}})^2 \to 0$ as $k \to \infty$.

{\bf Step 2:} We prove that the subsequence $(\mathcal{T}_{\ell_k})_{k \in \mathbb{N}_0}$ also guarantees that $\eta_{\ell_k}(\mathcal{M}_{\ell_k}) \to 0$ as $k \to \infty$. To this end, we first note that $T \in \mathcal{T}_{\ell_k} \backslash \mathcal{T}_{\ell_k+1}$ implies that $T$ is refined and hence $T \not \in \mathcal{T}_\infty$. Therefore, we see that
\begin{align*}
 \mathcal{M}_{\ell_k} \subseteq \mathcal{T}_{\ell_k} \backslash \mathcal{T}_{\ell_k+1}
 \subseteq \mathcal{T}_{\ell_k} \backslash \mathcal{T}_\infty
 = \mathcal{T}_{\ell_k} \backslash [\mathcal{T}_{\ell_k} \cap \mathcal{T}_\infty ]
 \reff{eq10:theoerem:estconv:new}= \mathcal{T}_{\ell_k} \backslash [\mathcal{T}_{\ell_{k+1}} \cap \mathcal{T}_{\ell_k}]
 = \mathcal{T}_{\ell_k} \backslash \mathcal{T}_{\ell_{k+1}}.
\end{align*}
This implies that $0 \le \eta_{\ell_k}(\mathcal{M}_{\ell_k}) \le \eta_{\ell_k}(\mathcal{T}_{\ell_k} \backslash \mathcal{T}_{\ell_{k+1}}) \to 0$
as $k \to \infty$.

{\bf Step 3:} For all fixed $\ell' \in \mathbb{N}_0$, we prove that
\begin{align}\label{eq:lemma70:new}
 \eta_{\ell_k}(\mathcal{T}_{\ell'} \cap \mathcal{T}_\infty) \to 0
 \quad \text{as } \ell' \le \ell_k \to \infty \text{ \  together with \ } k \to \infty.
\end{align}
To see this, we exploit the marking strategy~\eqref{eq:marking} and note with Step~2 that
\begin{align*}
 \max_{T \in \mathcal{T}_{\ell_k} \backslash \mathcal{M}_{\ell_k}} \eta_{\ell_k}(T)
 \le M \big( \eta_{\ell_k}(\mathcal{M}_{\ell_k}) \big) \to 0
 \quad \text{as } k \to \infty.
\end{align*}
For $\ell' \le \ell_k$, it holds that $\mathcal{T}_{\ell'} \cap \mathcal{T}_\infty \subseteq \mathcal{T}_{\ell_k} \backslash \mathcal{M}_{\ell_k}$ and hence
\begin{align*}
 \eta_{\ell_k}(T) \to 0 \quad \text{as } k \to \infty
 \quad \text{for all } T \in \mathcal{T}_{\ell'} \cap \mathcal{T}_\infty.
\end{align*}
Since $\mathcal{T}_{\ell'} \cap \mathcal{T}_\infty \subseteq \mathcal{T}_{\ell_k}$ is a fixed finite set, we conclude the proof of~\eqref{eq:lemma70:new}.

{\bf Step 4:} We prove that $(\mathcal{T}_{\ell_k})_{k \in \mathbb{N}_0}$ admits a subsequence $(\mathcal{T}_{\ell_{k_j}})_{j \in \mathbb{N}_0}$ such that $\eta_{\ell_{k_j}} \to 0$ as $j \to \infty$. To this end, let first $\mathcal{T}_H \in \mathbb{T}$ and $\mathcal{T}_h \in \mathbb{T}(\mathcal{T}_H)$.
Reduction~\eqref{eq:axiom:reduction}
proves that
\begin{subequations}\label{theoreom:new:eq1}
\begin{align}\nonumber
 \eta_h^2
 &=
 \eta_h(\mathcal{T}_h \backslash \mathcal{T}_H)^2 + \eta_h(\mathcal{T}_h \cap \mathcal{T}_H)^2
 \\&
 \le
 \contract{red} \, \eta_H(\mathcal{T}_H \backslash \mathcal{T}_h)^2 + \eta_h(\mathcal{T}_h \cap \mathcal{T}_H)^2
 + R(\norm{u_h - u_H}{\mathcal{X}})
 \\& \nonumber
 \le \contract{red} \, \eta_H^2 + \eta_h(\mathcal{T}_h \cap \mathcal{T}_H)^2
 + R(\norm{u_h - u_H}{\mathcal{X}}).
\intertext{If $\eta_H = 0$, stability~\eqref{eq:axiom:stability} and reduction~\eqref{eq:axiom:reduction} prove that}
 \eta_h^2 \nonumber
 &=
 \eta_h(\mathcal{T}_h \backslash \mathcal{T}_H)^2 + \eta_h(\mathcal{T}_h \cap \mathcal{T}_H)^2
 \\&
 \le \contract{red} \, \eta_H(\mathcal{T}_H \backslash \mathcal{T}_h)^2
 + 2 \, \eta_H(\mathcal{T}_H \cap \mathcal{T}_h)^2
 + [R + 2S^2](\norm{u_h - u_H}{\mathcal{X}})
 \\& \nonumber
 =  [R + 2S^2](\norm{u_h - u_H}{\mathcal{X}}),
\end{align}
\end{subequations}
where $[R+ 2 S^2](t) = R(t) + 2 S(t)^2$.
For the subsequence~$(\mathcal{T}_{\ell_k})_{k \in \mathbb{N}_0}$, we recall from~\eqref{eq10:theoerem:estconv:new} that $\mathcal{T}_{\ell_{k+n}} \cap \mathcal{T}_{\ell_k} = \mathcal{T}_{\ell_k} \cap \mathcal{T}_\infty$. Hence,
the estimates~\eqref{theoreom:new:eq1} read, for all $k\in \mathbb{N}_0$ and $n \in \mathbb{N}$,
\begin{align*}
 \eta_{\ell_{k+n}}^2
 &\le  \contract{red} \, \eta_{\ell_k}^2 + \eta_{\ell_{k+n}}(\mathcal{T}_{\ell_k} \cap \mathcal{T}_\infty)^2
 + R(\norm{u_{\ell_{k+n}} - u_{\ell_k}}{\mathcal{X}}),
 \quad \text{if $\eta_{\ell_k}^2 \neq 0$},
\intertext{resp.}
 \eta_{\ell_{k+n}}^2
 &\le [R+2S^2](\norm{u_{\ell_{k+n}} - u_{\ell_k}}{\mathcal{X}}),
 \quad \text{if $\eta_{\ell_k}^2 = 0$}.
\end{align*}
Let $0 < \contract{red} < \contract{red}' < 1$. Given $k \in \mathbb{N}_0$ with $\eta_{\ell_k}^2 \neq 0$, the convergence~\eqref{eq:lemma70:new}
allows us to pick some $n(k) \in \mathbb{N}$ such that
\begin{align*}
 \contract{red} \, \eta_{\ell_k}^2 + \eta_{\ell_{k+n(k)}}(\mathcal{T}_{\ell_k} \cap \mathcal{T}_\infty)^2
 \le  \contract{red}' \, \eta_{\ell_k}^2.
\end{align*}
In particular, we can choose a further subsequence $(\mathcal{T}_{\ell_{k_j}})_{j \in \mathbb{N}_0}$ of $(\mathcal{T}_{\ell_k})_{k \in \mathbb{N}_0}$ such that
\begin{align*}
 \eta_{\ell_{k_{j+1}}}^2
 \le \contract{red}' \, \eta_{\ell_{k_j}}^2 + [R+2S^2](\norm{u_{\ell_{k_{j+1}}} - u_{\ell_{k_j}}}{\mathcal{X}})
 \quad \text{for all } j \in \mathbb{N}_0.
\end{align*}
With $a_j = \eta_{\ell_{k_j}}^2$, $\rho = \contract{red}'$, and $b_j = [R+2S^2](\norm{u_{\ell_{k_{j+1}}} - u_{\ell_{k_j}}}{\mathcal{X}})$, {\sl a~priori} convergence~\eqref{assumption:limit_u} proves that
\begin{align*}
 0 \le a_{j+1} \le \rho \, a_j + b_j
 \quad \text{for all } j \in \mathbb{N}_0
 \quad \text{with } \lim_{j \to \infty} b_j = 0.
\end{align*}
By use of Lemma~\ref{axioms:lemma:cor:estimator_convergence}, we see that $a_j = \eta_{\ell_{k_j}}^2 \to 0$ as $j \to \infty$.

{\bf Step~5:} We prove that convergence of the subsequence $\eta_{\ell_{k_j}} \to 0$ already implies convergence of the full sequence $\eta_\ell \to 0$ as $\ell \to \infty$. To this end, we argue as in~\eqref{theoreom:new:eq1} and use $\contract{red} \le 2$ to see that, for all $\mathcal{T}_H \in \mathbb{T} \text{ and all } \mathcal{T}_h \in \mathbb{T}(\mathcal{T}_H)$,
\begin{align}\label{theoreom:new:eq2}
 \eta_h^2
 \le 2 \, \eta_H^2 + [R+2S^2] (\norm{u_h - u_H}{\mathcal{X}}).
\end{align}
Given $\varepsilon > 0$, there exists an index $\ell_{k_j}$ such that $\eta_{\ell_{k_j}}\le \varepsilon$ and $[R+2S^2](\norm{u_\ell - u_{\ell_{k_j}}}{\mathcal{X}}) \le \varepsilon$ for all $\ell \ge \ell_{k_j}$. For $\ell \ge \ell_{k_j}$, estimate~\eqref{theoreom:new:eq2} thus proves that
\begin{align*}
 \eta_\ell^2 \le 2\varepsilon^2 + \varepsilon.
\end{align*}
This concludes the proof.
\end{proof}

For the D\"orfler marking criterion~\eqref{eq:marking:doerfler}, the
refinement assumption~\eqref{eq:union_of_sons} exploited in the proof of Theorem~\ref{theoerem:estconv:new} can even be dropped. The following result (together with its very simple proof) is essentially the key argument in~\cite{afp2012}.

\begin{theorem}\label{theoerem2:estconv:new}
Suppose the properties~\eqref{eq:axiom:stability}--\eqref{eq:axiom:reduction} of the estimator. Consider the output of Algorithm~\ref{alg:abstract_algorithm}, where the marked elements $\mathcal{M}_\ell \subseteq \mathcal{T}_\ell$ satisfy the D\"orfler criterion~\eqref{eq:marking:doerfler} for some fixed marking parameter $0 < \theta \le 1$. Then, {\sl a~priori} convergence~\eqref{assumption:limit_u} yields estimator convergence~\eqref{eq:theoerem:estconv:new}.
\end{theorem}

\begin{proof}
Let $\mathcal{T}_H \in \mathbb{T}$ and $\mathcal{T}_h \in \mathbb{T}(\mathcal{T}_H)$. Arguing as for~\eqref{theoreom:new:eq1} and exploiting the Young inequality for some arbitrary $\delta > 0$ (instead of $\delta = 1$ above), we see that
\begin{align*}
 &\eta_h^2
 \le
 \contract{red} \, \eta_H(\mathcal{T}_H \backslash \mathcal{T}_h)^2 + (1+\delta) \, \eta_H(\mathcal{T}_h \cap \mathcal{T}_H)^2

 + [R + (1+\delta^{-1}) S^2] (\norm{u_h - u_H}{\mathcal{X}})
 \\& \quad
 = (1+\delta) \, \eta_H^2 - [(1+\delta) -\contract{red}] \, \eta_H(\mathcal{T}_H \backslash \mathcal{T}_h)^2

 + [R + (1+\delta^{-1}) S^2] (\norm{u_h - u_H}{\mathcal{X}}).
\end{align*}
For the sequence $(\mathcal{T}_\ell)_{\ell \in \mathbb{N}_0}$, the D\"orfler criterion~\eqref{eq:marking:doerfler} and $\mathcal{M}_\ell\subseteq\mathcal{T}_\ell\setminus\mathcal{T}_{\ell+1}$ yield that
\begin{align*}
 \theta \eta_\ell^2 \le \eta_\ell(\mathcal{M}_\ell)^2 \le \eta_\ell(\mathcal{T}_\ell \backslash \mathcal{T}_{\ell+1})^2.
\end{align*}
Combining the latter estimates, we obtain that
\begin{align*}
 \eta_{\ell+1}^2
 &\,\,\le\,\, (1+\delta) \, \eta_\ell^2 - [(1+\delta) -\contract{red}] \, \eta_\ell(\mathcal{T}_\ell \backslash \mathcal{T}_{\ell+1})^2
 +[R + (1+\delta^{-1}) S^2] (\norm{u_{\ell+1} - u_\ell}{\mathcal{X}}) \\
 &\reff{eq:marking:doerfler}\le \big((1+\delta) - [(1+\delta) -\contract{red}] \, \theta \big) \, \eta_\ell^2
 + [R + (1+\delta^{-1}) S^2] (\norm{u_{\ell+1} - u_\ell}{\mathcal{X}}).
\end{align*}
We define $\rho = (1+\delta) - [(1+\delta) -\contract{red}] \, \theta > 0$ as well as  $a_\ell = \eta_\ell^2$ and $b_\ell = [R + (1+\delta^{-1}) S^2] (\norm{u_{\ell+1} - u_\ell}{\mathcal{X}})$. Choosing $\delta > 0$ sufficiently small, we observe that $0 < \rho < 1$ and that
{\sl a~priori} convergence~\eqref{assumption:limit_u} proves that
\begin{align*}
 0 \le a_{\ell + 1} \le \rho \, a_\ell + b_\ell
 \quad \text{for all } \ell \in \mathbb{N}_0
 \quad \text{with } \lim_{\ell \to \infty} b_\ell = 0.
\end{align*}
By use of Lemma~\ref{axioms:lemma:cor:estimator_convergence}, we conclude that $a_\ell = \eta_\ell^2 \to 0$ as $\ell\to \infty$.
\end{proof}

\begin{remark}
The proof of Theorem~\ref{theoerem2:estconv:new} shows that in case of the D\"orfler marking~\eqref{eq:marking:doerfler}, it is sufficient to have stability~\eqref{eq:axiom:stability} and reduction~\eqref{eq:axiom:reduction} for one-level refinements, i.e., \eqref{eq:axiom:stability}--\eqref{eq:axiom:reduction} are only required for all $\mathcal{T}_H \in \mathbb{T}$, all $\mathcal{M}_H \subseteq \mathcal{T}_H$, and $\mathcal{T}_h = {\rm refine}(\mathcal{T}_H, \mathcal{M}_H)$.
\end{remark}%

\section{Examples}
\label{section:examples}

\input{01_obstacle}

\input{02_fractional}

\input{03_weaksing}

\input{04_hypsing}

\input{05_fembem}

\bibliographystyle{alpha}


\section{Introduction}

Let $A$ be an $n\times n$ matrix with eigenvalues $a_{1},\dots,a_{n}$.
It is well known that
\begin{align*}
a_{1}+a_{2}+\dots+a_{n} & =\mathrm{tr}(A)=A_{11}+A_{22}+\dots+A_{nn}\\
a_{1}a_{2}\cdots a_{n} & =\det(A)=\sum_{\sigma\in S_{n}}(-1)^{\sigma}A_{1,\sigma(1)}A_{2,\sigma(2)}\cdots A_{n,\sigma(n)}
\end{align*}
By applying the fundamental theorem of symmetric polynomials to the
characteristic polynomial of $A$, we find that there must be such
an equation between any symmetric polynomial in the eigenvalues of
$A$ and \emph{some} polynomial in the entries of $A$. We give an
explicit formula for this polynomial in terms of determinants, which
generalises the equations above to any symmetric polynomial:
\begin{align*}
\sum_{i\in\mathbb{N}^{n}}p_{i}a_{1}^{i_{1}}\cdots a_{n}^{i_{n}} & =\sum_{i\in\mathbb{N}^{n}}p_{i}\det(A_{1}^{i_{1}},\dots,A_{n}^{i_{n}})
\end{align*}
At first, this may seem surprising, since eigenvalues are independent
of the choice of basis, whereas taking $k$-th columns of a matrix
is clearly basis dependent. Indeed, each term $p_{i}\det(A_{1}^{i_{1}},\dots,A_{n}^{i_{n}})$
is on its own not basis independent, only the whole sum is, and only
if the $p_{i}$ are coefficients of a symmetric polynomial (which
means $p_{(i_{\sigma(1)},\dots,i_{\sigma(n)})}=p_{(i_{1},\dots,i_{n})}$
for any permutation $\sigma$, or $p_{i\circ\sigma}=p_{i}$ in short,
and that only finitely many $p_{i}$ are nonzero). More generally,
\begin{align}
\sum_{i\in\mathbb{N}^{n}}p_{i}\det(A_{1}^{(i_{1})},\dots,A_{n}^{(i_{n})})\label{eq:detAj}
\end{align}
is basis independent for any family of matrices $A^{(j)}$, not necessarily
powers $A^{j}$ of a single matrix. That is, if we substitue $A^{(j)}\mapsto S^{-1}A^{(j)}S$
for some invertible matrix $S$, its value does not change. Furthermore,
if the $A^{(j)}$ commute, we have the identity
\begin{align*}
\sum_{i\in\mathbb{N}^{n}}p_{i}a_{1}^{(i_{1})}\cdots a_{n}^{(i_{n})} & =\sum_{i\in\mathbb{N}^{n}}p_{i}\det(A_{1}^{(i_{1})},\dots,A_{n}^{(i_{n})})
\end{align*}
where $a_{k}^{(j)}$ is the $k$-th eigenvalue of $A^{(j)}$.

Our strategy to prove this is to define a quantity that is basis independent
because its definition makes no reference to a basis, and then show
that it is equal to (\ref{eq:detAj}). Once we have shown that (\ref{eq:detAj})
is invariant under basis transformations, we pick a basis in which
the determinants become products of eigenvalues, which is possible
if the $A^{(j)}$ commute.

By applying this identity to particular families of matrices, we get
the fundamental theorem of symmetric polynomials as a corollary, and
we are able to deduce various equations between eigenvalues and determinants,
such as those in the abstract.

\section{Preliminaries}

The proof is based on multilinear antisymmetric functions {[}1{]}.
\begin{defn*}
A function $f:V^{n}\to\mathbb{R}$ is
\end{defn*}
\begin{itemize}
\item Multilinear if $f$ is linear in each argument, i.e. $v_{k}\mapsto f(v_{1},\dots,v_{k},\dots,v_{n})$
is linear, with $v_{i}\in V$.
\item Antisymmetric if applying a permutation $\sigma\in S_{n}$ to its
arguments multiplies it by the sign of the permutation, i.e. $f(v_{\sigma(1)},\dots,v_{\sigma(n)})=(-1)^{\sigma}f(v_{1},\dots,v_{n})$,
with $v_{i}\in V$.
\end{itemize}
Multilinear antisymmetric functions form a vector space.
\begin{defn*}
Let $\bigwedge^{n}V^{*}\subset V^{n}\to\mathbb{R}$ be the space of
multilinear antisymmetric functions.
\end{defn*}
The only property of $\bigwedge^{n}V^{*}$ we will need is that $\dim(\bigwedge^{n}V^{*})=1$
if $n=\dim(V)$. In fact, in general $\dim(\bigwedge^{n}V^{*})={\dim(V) \choose n}$
{[}1{]}. Since it is the basis of our main theorem, we will give a
proof here.
\begin{lem*}
$\dim(\bigwedge^{n}V^{*})=1$ if $\dim(V)=n$.
\end{lem*}
\begin{proof}
Let $f\in\bigwedge^{n}V^{*}$ and $v_{1},\dots,v_{n}\in V$. Take
a basis $e_{1},\dots,e_{n}$ for $V$. Writing the $v_{i}$ in terms
of the basis, we have coefficients $a_{ij}\in\mathbb{R}$ for $i,j\in\{1\dots n\}$
such that $v_{i}=\sum_{j}a_{ij}e_{j}$. The proof proceeds by expanding
$f(v_{1},\dots,v_{n})$ by multilinearity, and then using the fact
that $f$ is zero whenever there is a duplicate argument, because
applying a permutation that swaps those two positions gives $f(\dots,w,\dots,w,\dots)=-f(\dots,w,\dots,w,\dots)$,
which implies $f(\dots,w,\dots,w,\dots)=0$. This allows us to convert
the sum over any pattern of indexing into a sum over permutations.
\begin{align*}
f(v_{1},\dots,v_{n}) & =f(\sum_{i_{1}=1}^{n}a_{1i_{1}}e_{i_{1}},\dots,\sum_{i_{n}=1}^{n}a_{1i_{n}}e_{i_{n}})\\
 & =\sum_{i_{1}=1}^{n}\dots\sum_{i_{n}=1}^{n}a_{1i_{1}}\dots a_{ni_{n}}f(e_{i_{1}},\dots,e_{i_{n}})\\
 & =\sum_{i\in\{1\dots n\}\to\{1\dots n\}}a_{1i(1)}\dots a_{ni(n)}f(e_{i(1)},\dots,e_{i(n)})\\
 & =\sum_{\sigma\in S_{n}}a_{1\sigma(1)}\dots a_{n\sigma(n)}f(e_{\sigma(1)},\dots,e_{\sigma(n)})\\
 & =\sum_{\sigma\in S_{n}}(-1)^{\sigma}a_{1\sigma(1)}\dots a_{n\sigma(n)}f(e_{1},\dots,e_{e})\\
 & =\det(a)f(e_{1},\dots,e_{n})
\end{align*}
Therefore, the value $f\in\bigwedge^{n}V^{*}$ is entirely determined
by its value on $f(e_{1},\dots,e_{n})$, and conversely, any value
chosen for $f(e_{1},\dots,e_{n})$ determines a multilinear antisymmetric
function $f\in\bigwedge^{n}V^{*}$ via the equation above, so the
space is one dimensional.
\end{proof}

\section{Proof}

Let $\mathbb{I}$ be some index set. We say that $p_{i}\in\mathbb{R}$
for $i\in\mathbb{I}^{n}$ are \emph{symmetric coefficients} if only
finitely many $p_{i}$ are nonzero, and $p_{i\circ\sigma}=p_{i}$
for all $i\in\mathbb{I}^{n}$ and permutations $\sigma\in S_{n}$.
Let $\vec{A}^{(j)}:V\to V$ for $j\in\mathbb{I}$ be a family of linear
maps on a vector space $V$ of dimension $n$.\\

\textcolor{teal}{}%
\noindent{\fboxsep 10pt\fbox{\begin{minipage}[t]{1\columnwidth - 2\fboxsep - 2\fboxrule}%
\begin{example*}
\textcolor{black}{A polynomial $p(x_{1},\dots,x_{n})$ on $n$ variables
can be written as a sum of monomials
\begin{align*}
p(x_{1},\dots,x_{n}) & =\sum_{i\in\mathbb{N}^{n}}p_{i}x_{1}^{i_{i}}\cdots x_{n}^{i_{n}}
\end{align*}
 where $p_{i}=p_{(i_{1},i_{2},\dots,i_{n})}$ is the coefficient of
$x_{1}^{i_{1}}x_{2}^{i_{2}}\dots x_{n}^{i_{n}}$. If $p$ is a symmetric
polynomial, then $p_{i}$ are symmetric coefficients with index set
$\mathbb{I}=\mathbb{N}$. For instance, if
\begin{align*}
p(x_{1},x_{2},x_{3}) & =x_{1}+x_{2}+x_{3}=x_{1}^{1}x_{2}^{0}x_{3}^{0}+x_{1}^{0}x_{2}^{1}x_{3}^{0}+x_{1}^{0}x_{2}^{0}x_{3}^{1}
\end{align*}
then $p_{i}=1$ for $i=(1,0,0),(0,1,0),(0,0,1)$ and $p_{i}=0$ otherwise.
Note that $p_{i\circ\sigma}=p_{i}$ for all permutations $\sigma\in S_{3}$.}
\end{example*}
\end{minipage}}}\textcolor{teal}{}\\

We now define a linear map $p(\vec{A}):\bigwedge^{n}V^{*}\to\bigwedge^{n}V^{*}$
in terms of the matrices $\vec{A}^{(j)}$ and symmetric coefficients
$p_{i}$. The key observation is that a linear map from a one-dimensional
vector space to itself is just multiplication by a scalar, so we are
in effect defining a scalar here. We shall soon see that this scalar
is precisely \ref{eq:detAj}.
\begin{defn}
\label{def:pA}If $p_{i}$ are symmetric coefficients, we define $p(\vec{A}):\bigwedge^{n}V^{*}\to\bigwedge^{n}V^{*}$
by
\begin{align*}
p(\vec{A})f(v_{1},\dots,v_{n}) & =\sum_{i\in\mathbb{I}^{n}}p_{i}f(A^{(i_{1})}v_{1},\dots,A^{(i_{n})}v_{n})
\end{align*}

This does indeed preserve antisymmetry, by change of summation variable
$i=j\circ\sigma$, and using $p_{j\circ\sigma}=p_{j}$:
\begin{align*}
p(\vec{A})f(v_{\sigma(1)},\dots,v_{\sigma(n)}) & =\sum_{i\in\mathbb{I}^{n}}p_{i}f(A^{(i_{1})}v_{\sigma(1)},\dots,A^{(i_{n})}v_{\sigma(n)})\\
 & =\sum_{j\in\mathbb{I}^{n}}p_{j\circ\sigma}f(A^{(j_{\sigma(1)})}v_{\sigma(1)},\dots,A^{(j_{\sigma(n)})}v_{\sigma(n)})\\
 & =\sum_{j\in\mathbb{I}^{n}}p_{j}(-1)^{\sigma}f(A^{(j_{1})}v_{1},\dots,A^{(j_{n})}v_{n})\\
 & =(-1)^{\sigma}p(\vec{A})f(v_{1},\dots,v_{n})
\end{align*}
\\
\end{defn}

\textcolor{teal}{}%
\noindent{\fboxsep 10pt\fbox{\begin{minipage}[t]{1\columnwidth - 2\fboxsep - 2\fboxrule}%
\begin{example*}
\textcolor{black}{Given a single matrix $A\in\mathbb{R}^{3\times3}$
and index set $\mathbb{I}=\mathbb{N}$, we can pick $\vec{A}^{(j)}=A^{j}$
to be the powers of that matrix. For the symmetric coefficients $p_{i}$
from the running example, we have
\begin{align*}
p(\vec{A})f(v_{1},v_{2},v_{3}) & =f(A^{1}v_{1},A^{0}v_{2},A^{0}v_{3})+f(A^{0}v_{1},A^{1}v_{2},A^{0}v_{3})+f(A^{0}v_{1},A^{0}v_{2},A^{1}v_{3})\\
 & =f(Av_{1},Iv_{2},Iv_{3})+f(Iv_{1},Av_{2},Iv_{3})+f(Iv_{1},Iv_{2},Av_{3})\\
 & =f(Av_{1},v_{2},v_{3})+f(v_{1},Av_{2},v_{3})+f(v_{1},v_{2},Av_{3})
\end{align*}
}
\end{example*}
\end{minipage}}}\textcolor{teal}{}\\

We use the notation $[p(\vec{A})]\in\mathbb{R}$ for the scalar corresponding
to $p(\vec{A})\in End(\bigwedge^{n}V^{*})$, so that $p(\vec{A})f=[p(\vec{A})]f$
for all $f\in\bigwedge^{n}V^{*}$. Since $p(\vec{A})$ has been defined
in terms of the $A^{(j)}$, the scalar $[p(\vec{A})]\in\mathbb{R}$
is a function of the maps $A^{(j)}:V\to V$. Since we have not used
the choice of a basis for $V$ to define $p(\vec{A})$, the scalar
$[p(\vec{A})]$ is manifestly invariant under change of basis.

Pick a basis $B:\mathbb{R}^{n}\to V$ for $V$ (with basis vectors
$b_{i}=Be_{i}$, where $e_{i}\in\mathbb{R}^{n}$ is the standard basis).
We have matrix representations $M^{(j)}=B^{-1}A^{(j)}B\in\mathbb{R}^{n\times n}$
for the $A^{(j)}:V\to V$, and we wish to calculate $[p(\vec{A})]\in\mathbb{R}$
explicitly terms of the entries of $M^{(j)}\in\mathbb{R}^{n\times n}$.
\begin{thm}
\label{thm:det-formula}$[p(\vec{A})]=\sum_{i\in\mathbb{I}^{n}}p_{i}\det(M_{1}^{(i_{1})},\dots,M_{n}^{(i_{n})})$\\
\end{thm}

\textcolor{teal}{}%
\noindent{\fboxsep 10pt\fbox{\begin{minipage}[t]{1\columnwidth - 2\fboxsep - 2\fboxrule}%
\begin{example*}
\textcolor{black}{For the running example, and $B=I$,
\begin{align*}
[p(\vec{A})] & =\det(A_{1},I_{2},I_{3})+\det(I_{1},A_{2},I_{3})+\det(I_{1},I_{2},A_{3})=A_{11}+A_{22}+A_{33}
\end{align*}
}
\end{example*}
\end{minipage}}}\textcolor{teal}{}\\

\begin{proof}
Since $p(\vec{A})$ is multiplication by a scalar $[p(\vec{A})]$,
\begin{align*}
p(\vec{A})f(v_{1},\dots,v_{n}) & =[p(\vec{A})]\cdot f(v_{1},\dots,v_{n})
\end{align*}
for all $f\in\bigwedge^{n}V^{*}$ and vectors $(v_{1},\dots,v_{n})$.
Taking $f(w_{1},\dots,w_{n})=\det(B^{-1}w_{1},\dots,B^{-1}w_{n})$
and $(v_{1},\dots,v_{n})=(Be_{1},\dots,Be_{n})$ to be the basis vectors,
on the right hand side
\begin{align*}
f(v_{1},\dots,v_{n}) & =\det(B^{-1}Be_{1},\dots,B^{-1}Be_{n})=\det(e_{1},\dots,e_{n})=1
\end{align*}
and on the left hand side
\begin{align*}
p(\vec{A})f(v_{1},\dots,v_{n}) & =\sum_{i\in\mathbb{I}^{n}}p_{i}\det(B^{-1}A^{(i_{1})}Be_{1},\dots,B^{-1}A^{(i_{n})}Be_{n})=\sum_{i\in\mathbb{I}^{n}}p_{i}\det(M_{1}^{(i_{1})},\dots,M_{n}^{(i_{n})})
\end{align*}
giving $[p(\vec{A})]=\sum_{i\in\mathbb{I}^{n}}p_{i}\det(M_{1}^{(i_{1})},\dots,M_{n}^{(i_{n})})$.
\end{proof}
This shows that the value of $[p(\vec{A})]=\sum_{i\in\mathbb{I}^{n}}p_{i}\det(M_{1}^{(i_{1})},\dots,M_{n}^{(i_{n})})$
does not depend on the basis, since the left hand side is defined
without reference to the basis. By picking a basis in which the $A^{(j)}$
are all upper triangular, we can relate $[p(\vec{A})]$ to the eigenvalues
of the $A^{(j)}$.
\begin{thm}
\label{thm:eigs-formula}If the $A^{(j)}$ commute, then $[p(\vec{A})]=\sum_{i\in\mathbb{I}^{n}}p_{i}a_{1}^{(i_{1})}\cdots a_{n}^{(i_{n})}$,
where $a_{k}^{(j)}$ are the eigenvalues of $A^{(j)}$.\\
\\
\textcolor{teal}{}%
\noindent{\fboxsep 10pt\fbox{\begin{minipage}[t]{1\columnwidth - 2\fboxsep - 2\fboxrule}%
\begin{example*}
\textcolor{black}{For the running example, suppose we have a basis
in which $A$ is upper triangular (e.g. the Jordan basis), then
\begin{align*}
[p(\vec{A})] & =\det(A_{1},I_{2},I_{3})+\det(I_{1},A_{2},I_{3})+\det(I_{1},I_{2},A_{3})=a_{1}+a_{2}+a_{3}
\end{align*}
 where $a_{1},a_{2},a_{3}$ are the eigenvalues of $A$.}
\end{example*}
\end{minipage}}}\textcolor{teal}{}\\
\end{thm}

\begin{proof}
Commuting matrices have a basis in which they are simultaneously upper
triangular, by Schur decomposition {[}2{]}. The diagonal of those
upper triangular matrices $M^{(j)}$ will contain the eigenvalues
$a_{k}^{(j)}$. In this case, $\det(M_{1}^{(i_{1})},\dots,M_{n}^{(i_{n})})$
is the determinant of an upper triangular matrix, with eigenvalues
$a_{k}^{(i_{k})}$ on the diagonal. Hence $\det(M_{1}^{(i_{1})},\dots,M_{n}^{(i_{n})})=a_{1}^{(i_{1})}\cdots a_{n}^{(i_{n})}$,
and substituting this into theorem (\ref{thm:det-formula}) gives
the eigenvalue formula for $p(\vec{A})$.
\end{proof}
The condition that the $A^{(j)}$ commute is not a necessary condition,
because there are upper triangular matrices that do not commute.

\section{Corollaries}

By combining theorems (\ref{thm:det-formula}) and (\ref{thm:eigs-formula})
we can justify the magic formula for converting symmetric expressions
involving eigenvalues into symmetric expressions involving determinants.
We first show this by example: let $A,B$ be commuting $3\times3$
matrices with eigenvalues $a_{1},a_{2},a_{3}$ and $b_{1},b_{2},b_{3}$.
Formally, we take the index set $\mathbb{I}=\{A,B\}$ in the theorems.
Next, we choose symmetric coefficients $p_{i}$ for $i\in\mathbb{I}^{3}$;
we must choose values for $p_{AAA},p_{AAB},p_{ABA},\dots,p_{BBB}$
that are symmetric under permutations of the indices. We choose $p_{AAB}=p_{ABA}=p_{BAA}=1$
and the rest $0$. The theorems (\ref{thm:det-formula}) and (\ref{thm:eigs-formula})
give us the equation
\begin{align*}
a_{1}a_{2}b_{3}+a_{1}b_{2}a_{3}+b_{1}a_{2}a_{3} & =\det(A_{1},A_{2},B_{3})+\det(A_{1},B_{2},A_{3})+\det(B_{1},A_{2},A_{3})
\end{align*}
We can now proceed to pick the matrices $A,B$ in this equation, as
long as they commute. For instance, given a matrix $C$, we could
pick $A=C^{3}$ and $B=C^{-1}$. Or, given commuting matrices $C,D$,
we could pick $A=C^{2}+D$ and $B=CD$. Or we could pick $A=(C+D)^{-1}$
and $B=\exp(CD)$. We thus obtain various relations between eigenvalues
and determinants by substituting, e.g. $A=(C+D)^{-1}$, $B=\exp(CD)$,
and $a_{i}=(c_{i}+d_{i})^{-1}$, $b_{i}=\exp(c_{i}d_{i})$ into the
equation.

Rather than trying to capture this general method in a theorem, we
present a few special cases.
\begin{cor}
Let $q(b_{1},\dots,b_{n})=\sum_{i\in\mathbb{N}^{n}}p_{i}b_{1}^{i_{1}}\dots b_{n}^{i_{n}}$
be a symmetric polynomial in the eigenvalues of an $n\times n$ matrix
$B$, then $q(b_{1},\dots,b_{n})=\sum_{i\in\mathbb{N}^{n}}p_{i}\det(B_{1}^{i_{1}},\dots,B_{n}^{i_{n}})$.
\end{cor}

\begin{proof}
Apply theorems (\ref{thm:det-formula}) and (\ref{thm:eigs-formula})
with index set $\mathbb{I}=\mathbb{N}$ and matrices $A^{(j)}=B^{j}$.
\end{proof}
\begin{cor}
Let $q(x)=\sum_{k=0}^{n}a_{k}x^{k}$ be a polynomial with roots $r_{k}$.
Then a symmetric polynomial in the roots $r_{k}$ can be written as
a polynomial in the coefficients $a_{k}$.
\end{cor}

\begin{proof}
Apply the previous corollary with $B$ being the companion matrix
of $q$. The eigenvalues of $B$ are the $r_{k}$. The entries of
the companion matrix are all $0$, or $1$, or $a_{k}$, so $\det(B_{1}^{i_{1}},\dots,B_{n}^{i_{n}})$
is a polynomial in the $a_{k}$.
\end{proof}
A symmetric polynomial $q$ can be seen as a function of a vector
$x\in\mathbb{R}^{n}$ satisfying $q(Px)=q(x)$ for permutation matrices
$P$.
\begin{defn}
A multivariate symmetric polynomial $q(X)$ is a polynomial function
of the entries of a matrix $X\in\mathbb{R}^{n\times m}$ satisfying
$q(PX)=q(X)$ for permutation matrices $P$.
\end{defn}

The permutation $P$ permutes the rows of $X$, but keeps each row
together. A symmetric polynomial is the special case $m=1$, when
each row consists of a single entry.
\begin{cor}
Let $B^{(j)}$ for $j\in\{1,\dots,m\}$ be commuting matrices with
eigenvalues $b_{i}^{(j)}$, and define $X_{ij}=b_{i}^{(j)}$ to be
the matrix of eigenvalues. Then a multivariate symmetric polynomial
$q(X)$ can be written as a polynomial in the entries of the $B^{(j)}$.
\end{cor}

\begin{proof}
Take $\mathbb{I}=\mathbb{N}^{m}$ and $A^{(g)}=(B^{(1)})^{g_{1}}\cdots(B^{(m)})^{g_{m}}$
in theorems (\ref{thm:det-formula}) and (\ref{thm:eigs-formula}).
\end{proof}

\section{Footnote}

The same definition (\ref{def:pA}) works for $p(\vec{A}):\bigwedge^{k}W^{*}\to\bigwedge^{k}V^{*}$
also when $W\neq V$ and $k\neq n$. We can view $p(\vec{A})f$ as
a generalised pullback of $f$ along a list of maps $A^{(j)}:V\to W$.
The ordinary pullback $A^{*}:\bigwedge^{k}W^{*}\to\bigwedge^{k}V^{*}$
is the special case of a single map $A:V\to W$. We have, in general
$p(X\vec{A}Y)=Y^{*}p(\vec{A})X^{*}$, where $X\vec{A}Y$ is simultaneous
conjugation $(X\vec{A}Y)^{(j)}=XA^{(j)}Y$. This gives us a slightly
stronger version of theorem (\ref{thm:det-formula}): if $k=\dim(V)=\dim(W)$,
then $X^{*}=\det(X)$ and $Y^{*}=\det(Y)$, so we see that if we multiply
the $M^{(j)}$ on the left by $X$ and on the right by $Y$, the value
of $[p(\vec{A})]$ gets multiplied by $\det(X)\det(Y)$. Theorem (\ref{thm:det-formula})
tells us that the value does not change if we do a basis transformation,
but here we get information about the case $X\neq Y^{-1}$. It is
also possible to generalise the proof of theorem (\ref{thm:det-formula})
directly.

\section*{References}

{[}1{]} Nicolas Bourbaki. Elements of mathematics, Algebra I. Springer-Verlag,
1989.

{[}2{]} R.A. Horn and C.R. Johnson. Matrix Analysis. Cambridge University
Press, 1985.
\end{document}