On some coupled local and nonlocal diffusion models

TL;DR

The paper develops energy-based models to couple local and nonlocal diffusion on a domain split by an interface, using a common weighted $H^1$ framework for the local part and weighted $H^s$ forms for the nonlocal part with $s\in(0,1)$. It introduces two energies, $E_I$ and $E_{II}$, and derives their Euler–Lagrange equations, revealing how transmission across the interface is governed by different nonlocal flux mechanisms depending on the chosen energy. A thorough regularity analysis is performed, particularly for $E_{II}$, and finite element discretizations are analyzed to obtain convergence rates under suitable Besov/Sobolev regularity; the role of the interface $\Sigma$ and graded meshes is emphasized. Numerical experiments illustrate interface flux behavior, coupling between disconnected subdomains, and the limiting behavior as $s\to1$, which recovers standard local transmission conditions, demonstrating the practical viability and theoretical depth of the coupled models.

Abstract

We study problems in which a local model is coupled with a nonlocal one. We propose two energies: both of them are based on the same classical weighted $H^1$-semi norm to model the local part, while two different weighted $H^s$-semi norms, with $s \in (0,1)$, are used to model the nonlocal part. The corresponding strong formulations are derived. In doing so, one needs to develop some technical tools, such as suitable integration by parts formulas for operators with variable diffusivity, and one also needs to study the mapping properties of the Neumann operators that arise. In contrast to problems coupling purely local models, in which one requires transmission conditions on the interface between the subdomains, the presence of a nonlocal operator may give rise to nonlocal fluxes. These nonlocal fluxes may enter the problem as a source term, thereby changing its structure. Finally, we focus on a specific problem, that we consider most relevant, and study regularity of solutions and finite element discretizations. We provide numerical experiments to illustrate the most salient features of the models.

Theorems & Definitions (1)

• Remark 1.1: issues when $s \in (0,1/2{}$] While in principle the energy $E_{I}$ is meaningful for $s\in (0,1)$, for $s \in (0,1/2]$ it has two issues: it does not admit a unique minimizer and it gives rise to two decoupled problems. For the sake of simplicity let us assume ${\sigma_{n\ell}} \equiv 1$ and $\sl \equiv 1$; in such a case, the energy above is $E_{I}[u] = \frac{1}{2} |u|^2_{H^1(\Omega_{1})}+ \frac{1}{2} |u|^2_{H^s(\Omega_{2})} - \int_\Omega f u \, dx.$ Regarding the first issue, because for $s \in (0,1/2]$ we have $H^s(\Omega_{2}) \equiv H^s_0(\Omega_{2})$, we observe that the energy-minimization problem would not be well-posed because solutions would be defined up to a constant on $\Omega_{2}$. Another way to see this is to notice that $u \mapsto \frac{1}{2} |u|^2_{H^s(\Omega_{2})}$ does not induce a norm in $H^s_0(\Omega_{2})$. The second issue is related to the fact that, for $s \in (0,1/2]$, the energy $E_I$ does not enforce any continuity on $\Sigma$. Therefore, when seeking for critical points of $E_I$ one can on the one hand minimize the functional $u \mapsto \frac{1}{2} |u|^2_{H^1(\Omega_{1})} - \int_{\Omega_{1}} f u \, dx$ over functions that vanish on $\partial \Omega$ and on the other hand seek critical points of the functional $u \mapsto \frac{1}{2} |u|^2_{H^s(\Omega_{2})} - \int_{\Omega_{2}} f u \, dx$. The second energy we inspect is $E_{II}[u]:= \frac{1}{2} \int_{\Omega_{1}} \sl |\nabla u |^2 dx+ \frac{C(n,s)}{4} \iint_{Q_{\Omega_{2}}} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))^2}{|x-y|^{n+2s}} dy dx - \int_\Omega f u \, dx ,$ where $Q_{\Omega_{2}} = ({{\mathbb R}^n} \times {{\mathbb R}^n}) \setminus (\Omega_{2}^c \times \Omega_{2}^c)$, and we also take the completion of ${\mathcal{D}}(\Omega)$ with respect to the norm associated with the Euler-Lagrange equations. With respect to the previous energy, $E_{II}$ now incorporates the interactions between $\Omega_{2}$ and $\Omega_{2}^c$. This gives rise to meaningful minimization problems for all $s \in (0,1)$, and this energy does not suffer any of the drawbacks mentioned in Remark \ref{['rem:issues']}. The presence of the integral over $\Omega_{1}\times\Omega_{2}$ corresponds to a 'flux density' between $\Omega_{1}$ and $\Omega_{2}$ that enforces continuity across the interface and effectively couples the problems even for $s \in (0,1/2]$. Both energies $E_I$ and $E_{II}$ induce Hilbert space norms in $\mathcal{D}(\Omega)$. The outline of the paper is as follows. In the next Section, we review two different definitions of fractional Laplacians, the related integration by part formulas and nonlocal Neumann operators. In Section \ref{['sec_EulerLagrange']}, we build the Euler-Lagrange equations in strong form for both energies $E_{I}$ and $E_{II}$. In particular, there is a difference in the way the transmission condition acts on the solutions, depending on the chosen energy. We then study the a priori regularity of the solution, focusing on the second energy $E_{II}$. In a last part (Section \ref{['sec:experiments']}), we carry out the numerical analysis for both energies, and provide some numerical results that highlight some salient features of our models. Along the manuscript, we use the notation $a \lesssim b$ to denote that $a \le C b$ with a constant $C$ that is independent of both $a$ and $b$; when relevant, we specify dependence. In this section, we write suitable integration by parts formulas for the weighted operators appearing in the energies we proposed in §\ref{['sec:l-nl-coupling']}. We shall make use of fractional-order Sobolev spaces and Besov spaces. We employ the notation from BoNo23 and refer to that work for elementary properties of these spaces. Just in this subsection, we let $s \in (0,1)$ be given, instead of $s \in (1/2,1)$ elsewhere. This work is concerned with fractional Laplacians with order $s$. In the whole space ${{\mathbb R}^n}$, there is a clear way to define such an operator: it is the pseudo-differential operator with symbols $|\xi|^{2s}$. Among several equivalent characterizations one can provide Kwasnicki2017, one can set for sufficiently smooth $u \colon {{\mathbb R}^n} \to {\mathbb R}$, $(-\Delta)^s u (x) := C(n,s) \hbox{p.v.} \int_{{\mathbb R}^n} \frac{u(x)-u(y)}{|x-y|^{n+2s}} dy, \quad x \in {{\mathbb R}^n},$ with $C(n,s) := \frac{2^{2s} s \Gamma(s+\frac{n}{2})}{\pi^{n/2} \Gamma(1-s)}.$ The operator \ref{['eq:def-ifl']} arises as the infinitesimal generator of a $2s$-stable Lévy process, and therefore has been widely employed in a number of applications, including the modeling of market fluctuations ContTankov and subsurface flow through porous media BensonWheatcraft. Our interest in the fractional Laplacian and related operators stems from the analysis of interface problems involving dielectrics and metamaterials BoCi17, where the relaxation of the classical $H^1$-setting by introducing $H^s$-forms --and consequently leading to the use of the fractional Laplacian \ref{['eq:def-ifl']}-- yields a reduction of the so-called critical interval and therefore to a more stable behavior of solutions. There is not a unique way to define fractional Laplacians on arbitrary domains $\Lambda \subset {{\mathbb R}^n}$. We briefly comment on three non-equivalent definitions. The spectral fractional Laplacian corresponds to a non-integer power of the Laplace operator with Dirichlet boundary conditions in the spectral sense. A second approach consists in using definition \ref{['eq:def-ifl']}; we refer to such an operator as integral fractional Laplacian. As opposed to the spectral operator, the integral one does not depend on the domain $\Lambda$; this definition keeps the probabilistic interpretation of the fractional Laplacian defined over ${{\mathbb R}^n}$. However, one should also observe that \ref{['eq:def-ifl']} involves the values of $u$ outside of $\Lambda$. A third way to define a fractional Laplacian on a domain $\Lambda$ consists in restricting the integral in \ref{['eq:def-ifl']} to the domain $\Lambda$ itself. Namely, we define the censored fractional Laplacian (or regional fractional Laplacian) of order $s$ of a function $u \colon \Lambda \to {\mathbb R}$ by $(-\Delta)_\Lambda^s u (x) := C(n,s) \hbox{p.v.} \int_\Lambda \frac{u(x)-u(y)}{|x-y|^{n+2s}} dy, \quad x \in {{\mathbb R}^n}.$ This operator arises from taking first variation of the functional $u \mapsto \frac{1}{2}|u|^2_{H^s(\Lambda)}$. While in principle this operator is meaningful for $s \in (0,1)$, the case $s>1/2$ is the most interesting: in such a case, the Dirichlet problem for this operator involves boundary values on $\partial \Lambda$. This work is concerned with operators related to \ref{['eq:def-ifl']} and \ref{['eq:def-rfl']} on bounded, Lipschitz domains in ${{\mathbb R}^n}$. These operators shall arise in the Euler-Lagrange equations related to energies \ref{['eq:energy-I']} and \ref{['eq:energy-III']}, respectively. We let ${\sigma_{n\ell}}$ be as in the beginning of Section \ref{['sec:l-nl-coupling']}. By analogy with \ref{['eq:def-ifl']}, let us consider the weighted integral fractional Laplacian of order $s$, $(-\Delta_{\sigma_{n\ell}})^s u (x) := C(n,s) \hbox{p.v.}\int_{{{\mathbb R}^n}} {\sigma_{n\ell}}(x,y)\frac{u(x)-u(y)}{|x-y|^{n+2s}} dy.$ The symmetry of ${\sigma_{n\ell}}$ is required to have a problem with a variational structure; if one starts from an energy with non-symmetric diffusivity, then the associated operator in the Euler-Lagrange equations corresponds to a symmetrization of such a diffusivity. We can also consider a variant of \ref{['eq:def-rfl']} with variable diffusivity, the weighted censored fractional Laplacian of order $s$, $(-\Delta_{\sigma_{n\ell}})^s_\Lambda u (x) := C(n,s) \int_{\Lambda} {\sigma_{n\ell}}(x,y)\frac{u(x)-u(y)}{|x-y|^{n+2s}} dy, \quad x \in \Lambda.$ While operators like \ref{['eq:def-wifl']} and \ref{['eq:def-wrfl']} are well-defined under much more general assumptions on ${\sigma_{n\ell}}$, in practice we will require certain additional smoothness on ${\sigma_{n\ell}}$ in order to have integration by parts formulas or to exploit regularity estimates for the solutions of the related Dirichlet problems. We revisit this point in Section \ref{['sec:regularity']} below. Finally, we note that the pointwise definition of both \ref{['eq:def-wifl']} and \ref{['eq:def-wrfl']} is meaningful for sufficiently smooth functions. On the one hand, we have $(-\Delta_{\sigma_{n\ell}})^s : \mathcal{S}({{\mathbb R}^n}) \to \mathcal{S}({{\mathbb R}^n})$, and we can extend the definition to tempered distributions in a standard fashion, namely $(-\Delta_{\sigma_{n\ell}})^s : \mathcal{S}'({{\mathbb R}^n}) \to \mathcal{S}'({{\mathbb R}^n})$ by $\langle (-\Delta_{\sigma_{n\ell}})^s T , {\varphi} \rangle := \langle T , (-\Delta_{\sigma_{n\ell}})^s {\varphi} \rangle, \quad \forall T \in \mathcal{S}'({{\mathbb R}^n}), \, {\varphi} \in \mathcal{S}({{\mathbb R}^n}).$ On the other hand, the operator $(-\Delta_{\sigma_{n\ell}})^s_\Lambda$ is nonlocal and does not map neither $\mathcal{S}({{\mathbb R}^n})$ nor $\mathcal{D}(\Lambda)$ onto themselves. For the sake of our analysis, we simply point out that if ${\varphi} \in \mathcal{D}(\Lambda)$ then $(-\Delta_{\sigma_{n\ell}})^s_\Lambda {\varphi} \in L^2(\Lambda)$; additionally, integrating by parts (see \ref{['eq:ibp-censored-with-sigma']} below) we have $( (-\Delta_{\sigma_{n\ell}})^s_\Lambda u , {\varphi} ) = (u , (-\Delta_{\sigma_{n\ell}})^s_\Lambda {\varphi} ), \quad \forall u , {\varphi} \in \mathcal{D}(\Lambda).$ Therefore we can extend the definition of $(-\Delta_{\sigma_{n\ell}})^s_\Lambda \colon L^2(\Lambda) \to \mathcal{D}'(\Lambda)$ by $\langle (-\Delta_{\sigma_{n\ell}})^s_\Lambda u , {\varphi} \rangle := (u , (-\Delta_{\sigma_{n\ell}})^s_\Lambda {\varphi} ), \quad \forall u \in L^2(\Lambda), \, {\varphi} \in \mathcal{D}(\Lambda).$ We next review some integration by parts formulas for the (weighted) integral and censored fractional Laplacians on bounded domains. For the integral operator \ref{['eq:def-ifl']}, DiRoVa17 derives an identity of this type. A straightforward adaptation of the proof to the variable diffusivity setting gives rise to the following. Let $\Lambda \subset {\mathbb R}^n$ be bounded and Lipschitz, ${\sigma_{n\ell}}\in L^\infty(Q_\Lambda)$ be symmetric, and bounded functions $u,v \in C_c^2({\mathbb R}^n)$, $C(n,s)\iint_{Q_\Lambda} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx=\int_\Lambda v (-\Delta_{\sigma_{n\ell}})^s u\,dx + \int_{\Lambda^c} v\, \mathcal{N}_{s,{\sigma_{n\ell}}}^\Lambda u \,dx,$ where $\mathcal{N}_{s,{\sigma_{n\ell}}}^\Lambda u (x) := C(n,s) \int_\Lambda {\sigma_{n\ell}}(x,y) \frac{u(x)-u(y)}{|x-y|^{n+2s}} dy, \quad x \in \Lambda^c$ is a weighted nonlocal Neumann operator. The operator \ref{['eq:def-wNs']}, with ${\sigma_{n\ell}} \equiv 1$, has been proposed in DiRoVa17; while there is no widely accepted definition of a Neumann condition for the integral fractional Laplacian \ref{['eq:def-ifl']}, the one we chose has the advantage of giving rise to a variational structure while keeping the probabilistic interpretation of a Neumann condition as a random reflection --according to a Lévy flight-- of a particle inside the domain. We refer to DiRoVa17 for a thorough comparison between this and other existing approaches. As for the censored fractional Laplacian, according to Guan and GalWarma (see also Warma for a more general version), we have the following integration by parts formula. Given a bounded $C^{2}$ domain $\Lambda$, if $s \in (1/2,1)$, $v \in H^s(\Lambda)$, and $u \in C(\overline\Lambda)$ can be written as $u (x) = \phi(x) d(x,\partial\Lambda)^{2s-1} + \psi(x)$ for all $x\in \Lambda$, with $\phi,\psi \in C^2(\overline\Lambda)$, then $\frac{C(n,s)}{2}\iint_{\Lambda \times \Lambda} \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx=\int_\Lambda v (-\Delta)^s_\Lambda u \,dx + \int_{\partial\Lambda} v \mathcal{M}_s^\Lambda u \,d\sigma.$ Above, $\mathcal{M}_s^\Lambda$ is the fractional normal derivative $\mathcal{M}_s^\Lambda u (x) := C'(s) \lim_{t\to 0^+} \frac{u(x-t{\mathbf n})-u(x)}{t^{2s-1}}, \quad x \in \partial\Lambda,$ where ${\mathbf n}$ is the outward normal on $\partial\Lambda$, and $C'(s)$ is the constant introduced in Warma or Guan, $C'(s) := \frac{C(n,s)}{2s(2s-1)} \left(\int_0^\infty \frac{|t-1|^{1-2s} - \max\{1,t\}^{1-2s}}{t^{2-2s}} \, dt \right) \left(\int_{ \{|x|=1, x_n > 0\}} x_n^{2s} \, dx \right).$ We are interested in the extension of \ref{['eq:ibp-censored']} to the variable diffusivity case. Reference Guan provides sufficient conditions for such an integration by parts formula to hold. If ${\sigma_{n\ell}} \in C^1(\overline\Lambda)$, then taking $\psi_1 = {\sigma_{n\ell}}$, $\psi_2 = 0$ in that theorem allows one to get the formula below with $\mathcal{M}_{s,{\sigma_{n\ell}}}^{\Lambda} u (x) := {\sigma_{n\ell}}(x,x) \mathcal{M}_s^{\Lambda} u (x),$ namely $\frac{C(n,s)}{2}\iint_{\Lambda\times\Lambda} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx=\int_{\Lambda} v (-\Delta_{\sigma_{n\ell}})^s_{\Lambda} u \,dx+ \int_{\partial{\Lambda}} v \mathcal{M}_{s,{\sigma_{n\ell}}}^{\Lambda} u \,d\sigma.$ In view of the integration by parts formula \ref{['eq:wibp']}, to realize the mapping properties of the weighted Neumann operator $\mathcal{N}_{s,{\sigma_{n\ell}}}^\Lambda u$ it suffices to characterize the trace space associated to the energy. Such an effort has been carried out, for instance, in BogdanGrHe22 by using robust Poisson kernel estimates. Here, we follow reference GrKa23 and introduce the space $\mathcal{T}^s(\Lambda^c) := \{ g \colon \Lambda^c \to {\mathbb R} \hbox{measurable} \colon \| g \|_{\mathcal{T}^s(\Lambda^c)} < \infty \},$ where $\| g \|_{\mathcal{T}^s(\Lambda^c)}:= \left( \| g \|_{L^2(\Lambda^c; \mu_s)}^2 + | g |_{\mathcal{T}^s(\Lambda^c)}^2\right)^{1/2} ,| g |_{\mathcal{T}^s(\Lambda^c)}:= \left(\iint_{\Lambda^c \times \Lambda^c} \frac{(g(x) - g(y))^2}{((|x-y| + d_x + d_y ) \wedge 1)^n} \, \mu_s(dx) \mu_s(dy)\right)^{1/2},$ $d_x = \hbox{dist}(x, \partial \Lambda)$, and $\mu_s(dx) := \frac{(1-s) \chi_{\Lambda^c}(x)}{d_x^s (1+d_x)^{n+s}} \, dx.$ To deal with the natural space for the energy $E_{II}$, it shall be convenient to introduce the space $V^s(\Lambda | {\mathbb R}^n) := \{ u \colon {\mathbb R}^n \to {\mathbb R} \hbox{measurable} \colon \| u \|_{V^s(\Lambda | {\mathbb R}^n)} < \infty \},$ endowed with $\| u \|_{V^s(\Lambda | {\mathbb R}^n)}:= (\| u \|_{L^2(\Lambda)}^2 + |u|_{V^s(\Lambda | {\mathbb R}^n)}^2)^{1/2},|u|_{V^s(\Lambda | {\mathbb R}^n)}:= \left(\frac{C(n,s)}{2} \iint_{Q_\Lambda}\frac{(u(x)-u(y))^2}{|x-y|^{n+2s}} dy dx\right)^{1/2}.$ Let $t \in [s, \min\{1,2s\})$. Then, a straightforward application of the Cauchy-Schwarz inequality yields $\frac{C(n,s)}{2} \iint_{Q_\Lambda} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx \lesssim |u|_{V^t(\Lambda | {\mathbb R}^n)} |v|_{V^{2s-t}(\Lambda | {\mathbb R}^n)}$ for all $u \in V^{t}(\Lambda | {\mathbb R}^n), \, v \in V^{2s-t}(\Lambda | {\mathbb R}^n)$, with a constant depending on $n,s,t, {\sigma_{n\ell}}$. With the notation above, we have the following trace (restriction) and lifting (extension) result GrKa23. Let $\Lambda$ be a bounded Lipschitz domain and $t \in (0,1)$. There exists a continuous trace operator ${\rm Tr}_t \colon V^t(\Lambda | {\mathbb R}^n) \to \mathcal{T}^t(\Lambda^c)$.There exists a continuous lifting operator ${\rm Ext}_t \colon \mathcal{T}^t(\Lambda^c) \to V^t(\Lambda | {\mathbb R}^n)$. In both cases, the continuity constants depend on $\Omega$ and on a lower bound on $t$, namely, they are uniformly bounded as $t \to 1$. We can combine Theorem \ref{['thm:trace-extension']} with the integration by parts formula \ref{['eq:wibp']} to extend the operator $\mathcal{N}_{s,{\sigma_{n\ell}}}^\Lambda$ to a broader class of functions than $C^2_c({\mathbb R}^n)$. Indeed, let $t \in [s, \min\{1,2s\})$ and $\tau \in [0, \min\{1/2, 2s - t\} )$. We define the space $\mathcal{H}^t(\Lambda)$ through $\mathcal{H}^t(\Lambda) = \overline{C^\infty_c({\mathbb R}^n)}^{\| \cdot \|_{\mathcal{H}^t(\Lambda)}}, \quad \| u \|_{\mathcal{H}^t(\Lambda)} := \| u \|_{V^t(\Lambda | {\mathbb R}^n)} + \| (-\Delta_{\sigma_{n\ell}})^s u \|_{H^{-\tau}(\Lambda)}.$ Here, because $\tau<1/2$, we understand the restriction of the distribution $(-\Delta_{\sigma_{n\ell}})^s u$ in the dual space of $\widetilde{H}^\tau(\Lambda)$, namely, acting on functions that are extended by zero on $\Lambda^c$. We let $L \colon \mathcal{H}^t(\Lambda) \times V^{2s-t}(\Lambda | {\mathbb R}^n) \to {\mathbb R}$, $L(u,v) :=\frac{C(n,s)}{2} \iint_{Q_\Lambda} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx- \int_\Lambda v (-\Delta_{\sigma_{n\ell}})^s u\,dx .$ By Remark \ref{['rem:Cauchy-Schwarz']}, $L$ is a bounded, linear operator: $|L(u,v)|\lesssim \left( |u|_{V^t(\Lambda | {\mathbb R}^n)} |v|_{V^{2s-t}(\Lambda | {\mathbb R}^n)} + \| v \|_{H^\tau(\Lambda)}\| (-\Delta_{\sigma_{n\ell}})^s u \|_{H^{-\tau}(\Lambda)} \right)\lesssim \|u \|_{\mathcal{H}^t(\Lambda)} \|v\|_{V^{2s-t}(\Lambda | {\mathbb R}^n)},$ because $\| v \|_{H^\tau(\Lambda)} = (\| v \|_{L^2(\Lambda)}^2 + | v |_{H^\tau(\Lambda)}^2)^{1/2} \lesssim \|v\|_{V^{2s-t}(\Lambda | {\mathbb R}^n)}$ since $\tau \le 2s-t$ and $\Lambda$ is a bounded, Lipschitz set. Next, if $u \in C^2_c({\mathbb R}^n)$ and $v \in \mathcal{D}(\Lambda)$, we can use the integration by parts formula from Lemma \ref{['lemma:ibp-snl']} to observe $L(u,v) =\int_{\Lambda^c} v \mathcal{N}_s^\Lambda u\,dx = 0$ and, by density, we deduce that $L(u,v) = 0$ for all $u \in \mathcal{H}^t(\Lambda)$, $v \in \widetilde{H}^{2s-t}(\Lambda)$. We define $L_{\Lambda^c} \colon \mathcal{H}^t(\Lambda) \times \mathcal{T}^{2s-t}(\Lambda^c)$, $L_{\Lambda^c}(u,g) = L(u, {\rm Ext}_{2s-t}(g) ),$ where ${\rm Ext}_{2s-t} \colon \mathcal{T}^{2s-t}(\Lambda^c) \to V^{2s-t}(\Lambda | {\mathbb R}^n)$ is an extension operator as in Theorem \ref{['thm:trace-extension']}. It is clear that $L_{\Lambda^c}$ does not depend on how we extend $g$ to $\Lambda$: any other extension ${\rm Ext'}(g) \in V^{2s-t}(\Lambda | {\mathbb R}^n)$ would yield ${\rm Ext}_{2s-t}(g) - {\rm Ext'}(g) \in \widetilde{H}^{2s-t}(\Lambda)$ and therefore $L(u, {\rm Ext'}(g)) = L(u, {\rm Ext}_{2s-t}(g))$. Given $u \in \mathcal{H}^t(\Lambda)$, $g \in \mathcal{T}^{2s-t}(\Lambda^c)$, we then have by \ref{['eq:L-bounded']}, $|L_{\Lambda^c}(u,g)| \lesssim \| u \|_{\mathcal{H}^t(\Lambda)} \| {\rm Ext}_{2s-t}(g) \|_{V^{2s-t}(\Lambda | {\mathbb R}^n)} \lesssim \| u \|_{\mathcal{H}^t(\Lambda) } \| g \|_{ \mathcal{T}^{2s-t}(\Lambda^c)}.$ The preceding discussion leads us to defining, for any $t \in [s, \min \{1, 2s\})$, $\mathcal{N}_{s,{\sigma_{n\ell}}}^\Lambda \colon \mathcal{H}^t(\Lambda) \to \left( \mathcal{T}^{2s-t}(\Lambda^c)\right)'$ by $\langle \mathcal{N}_{s,{\sigma_{n\ell}}}^\Lambda u, g \rangle := L_{\Lambda^c}(u,g).$ Naturally, by construction this is a bounded operator; additionally, if $u$ is sufficiently smooth we can identify $\mathcal{N}_{s,{\sigma_{n\ell}}}^\Lambda$ with its pointwise definition, and the integration by parts formula remains valid. We introduced a parameter $t \in [s, \min \{1, 2s\})$ in the previous definition of the nonlocal Neumann operator $\mathcal{N}_{s,{\sigma_{n\ell}}}^\Lambda$. Unlike Sobolev spaces defined on bounded domains, the spaces $V^s(\Lambda|{\mathbb R}^n)$ are not ordered with respect to $s$: in general, it is not true that $V^t(\Lambda|{\mathbb R}^n) \subset V^s(\Lambda|{\mathbb R}^n)$ if $t > s$ because the former space can accommodate functions with a stronger growth at infinity than the latter. Therefore, formally, the operator \ref{['eq:def-Nsnl']} depends on $t$. However, for the sake of this work, if we restrict our attention to functions that have bounded support, then $t$ plays no significant role as all the $t$-dependent operators coincide with the pointwise definition \ref{['eq:def-wNs']} for functions in $C^\infty_c({\mathbb R}^n)$. To define the nonlocal Neumann operator $\mathcal{N}_{s,{\sigma_{n\ell}}}^\Lambda u$, we not only required the function $u$ to belong to the natural energy space $V^s(\Lambda | {\mathbb R}^n)$ but also that $(-\Delta_{\sigma_{n\ell}})^s u \in H^{-\tau}(\Lambda)$ for some $\tau < 1/2$. This is in alignment with the classical (local) case, in which one cannot define normal derivatives of functions on $H^1(\Omega)$ but requires some integrability of their Laplacian (see, for example, Sayas-book). For the two energies we introduced in §\ref{['sec:l-nl-coupling']}, that extend \ref{['eq:local_energy']} to a coupled local/nonlocal setting, we discuss the strong form of the resulting Euler-Lagrange equations. We recall that $\sl, {\sigma_{n\ell}}$ are bounded and uniformly positive functions. Conceptually, we follow a standard procedure: we first consider the first variation of the energies with respect to smooth functions supported in either $\Omega_{1}$, $\Omega_{2}$, to obtain suitable integro-differential PDEs on such domains, and afterwards with respect to smooth functions whose support intersects $\Sigma$ to derive transmission conditions. In doing so, we find some remarkable differences between these energies. We begin with the energy $E_{I}$ in \ref{['eq:energy-I']}. We assume $\Omega_{2}$ to be a $C^2$ domain, so that the integration by parts formula \ref{['eq:ibp-censored-with-sigma']} holds. The Euler-Lagrange equations associated to $E_I$ write: find $u\in \mathcal{H}_{I}$ such that $\int_{\Omega_{1}} \sl \nabla u \cdot \nabla v\, dx + \frac{ C(n,s) }{2} \iint_{\Omega_{2}^2} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx= \int_\Omega f v \, dx, \qquad \forall v \in \mathcal{H}_{I},$ where $\mathcal{H}_{I}$ is the completion of ${\mathcal{D}}(\Omega)$ with respect to the norm associated with the energy, $v \ \mapsto\ \int_{\Omega_{1}} \sl |\nabla v|^2\, dx+ \left( \frac{ C(n,s) }{2} \iint_{\Omega_{2}\times\Omega_{2}} {\sigma_{n\ell}}(x,y) \frac{(v(x)-v(y))^2}{|x-y|^{n+2s}} dy dx \right)^{1/2}.$ Step 1. We take $v\in{\mathcal{D}}(\Omega_{1})$ in \ref{['eq:Euler-Lagrange-E_I']} and integrate by parts to find that $\langle -{\rm div}\, (\sl \nabla u), v \rangle = \int_{\Omega_{1}} \sl \nabla u \cdot \nabla v\, dx = \int_{\Omega_{1}} f v \, dx \qquad \forall v \in {\mathcal{D}}(\Omega_{1}).$ Thus, in ${\mathcal{D}}'(\Omega_{1})$, and even in $L^2(\Omega_{1})$ since $f_{|\Omega_{1}}\in L^2(\Omega_{1})$, $-{\rm div}\, (\sl \nabla u) = f \hbox{in}\Omega_{1}.$ Step 2. Next, we take $v\in{\mathcal{D}}(\Omega_{2})$ in \ref{['eq:Euler-Lagrange-E_I']} to obtain, for all $v \in {\mathcal{D}}(\Omega_{2})$, $\frac{C(n,s)}{2} \iint_{\Omega_{2}\times\Omega_{2}} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx = \int_{\Omega_{2}} f v \, dx.$ On the other hand, using \ref{['eq:def-delta_domain-dist']} and the fact that $\mathcal{H}_{I} \subset H^s(\Omega_{2})$ so that we can apply \ref{['eq:ibp-censored-with-sigma']}, we have $\frac{ C(n,s) }{2} \iint_{\Omega_{2}\times\Omega_{2}} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx = \langle (-\Delta_{\sigma_{n\ell}})^s_{\Omega_{2}} u, v \rangle.$ Thus, in ${\mathcal{D}}'(\Omega_{2})$, and even in $L^2(\Omega_{2})$ since $f_{|\Omega_{2}}\in L^2(\Omega_{2})$, $(-\Delta_{\sigma_{n\ell}})^s_{\Omega_{2}} u = f \hbox{in}\Omega_{2}.$ Step 3. Finally, let us take $v\in{\mathcal{D}}(\Omega)$ in \ref{['eq:Euler-Lagrange-E_I']}. Formula \ref{['eq:u-Hdiv']} yields $\sl \nabla u |_{\Omega_{1}} \in \mathbf{H}({\rm div}\,,\Omega_{1})$, and we therefore have $\int_{\Omega_{1}} \sl \nabla u \cdot \nabla v\, dx= - \int_{\Omega_{1}} {\rm div}\,(\sl \nabla u ) v\, dx + \langle (\sl \frac{\partial u}{\partial {\mathbf n}})_{|\Omega_{1}}, v \rangle_{\widetilde{H}^{1/2}(\Sigma)}= \int_{\Omega_{1}} f v\, dx + \langle (\sl \frac{\partial u}{\partial {\mathbf n}})_{|\Omega_{1}}, v \rangle_{\widetilde{H}^{1/2}(\Sigma)},$ where $\widetilde{H}^{1/2}(\Sigma):= \{ g \in H^{1/2}(\partial \Omega_{1}) \colon g |_{{\Gamma_1}} = 0 \} = \{ g \in H^{1/2}(\partial \Omega_{2}) \colon g |_{{\Gamma_2}} = 0 \}$, and the equality holds algebraically and topologically. Since $(-\Delta_{\sigma_{n\ell}})^s_{\Omega_{2}} u\in L^2(\Omega_2)$, we can apply the integration by parts formula \ref{['eq:ibp-censored-with-sigma']} and replace in \ref{['eq:Euler-Lagrange-E_I']} to arrive at $\langle (\sl \frac{\partial u}{\partial {\mathbf n}})_{|\Omega_{1}}, v \rangle_{\widetilde{H}^{1/2}(\Sigma)} + \int_{\Sigma} v \mathcal{M}_{s,{\sigma_{n\ell}}}^{\Omega_2} u \,d\sigma = 0, \quad \forall v \in {\mathcal{D}}(\Omega) .$ In other words, we have shown that $(\sl \frac{\partial u}{\partial {\mathbf n}})_{|\Omega_{1}} + \mathcal{M}_{s,{\sigma_{n\ell}}}^{\Omega_2} u = 0 \hbox{on} \Sigma.$ To conclude, the minimization problem associated with the energy $E_I$ corresponds to the following: find $u\in \mathcal{H}_{I}$ such that $\left\lbrace -{\rm div}\, (\sl \nabla u) = f\hbox{in} \Omega_{1},(-\Delta_{\sigma_{n\ell}})^s_{\Omega_{2}} u = f\hbox{in} \Omega_{2},(\sl \frac{\partial u}{\partial {\mathbf n}})_{|\Omega_{1}} + (\mathcal{M}_{s,{\sigma_{n\ell}}}^{\Omega_2} u)_{|\Omega_{2}} = 0\hbox{on} \Sigma. \right.$ Thus, for this energy the coupling between the local and the nonlocal models occurs through the interface $\Sigma$. Remarkably, the resulting transmission condition involves the fractional normal derivative operator $\mathcal{M}_{s,{\sigma_{n\ell}}}^{\Omega_2}$. For the energy $E_{II}$ in \ref{['eq:energy-III']}, the Euler-Lagrange equations write: find $u\in \mathcal{H}_{II}$ such that $\int_{\Omega_{1}} \sl \nabla u \cdot \nabla v\, dx + \frac{C(n,s)}{2} \iint_{Q_{\Omega_{2}}} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx= \int_\Omega f v \, dx, \quad \forall v \in \mathcal{H}_{II},$ where $\mathcal{H}_{II}$ is the completion of ${\mathcal{D}}(\Omega)$ with respect to the norm associated with the Euler-Lagrange equations $v \ \mapsto\ \left( \int_{\Omega_{1}} \sl |\nabla v|^2\, dx + \frac{C(n,s)}{2} \iint_{Q_{\Omega_{2}}} {\sigma_{n\ell}}(x,y) \frac{(v(x)-v(y))^2}{|x-y|^{n+2s}} dy dx\right)^{1/2}$ which, since $\sigma_{\max}\ge\sl, {\sigma_{n\ell}} \ge\sigma_{\min}>0$ a.e. and using the notation \ref{['eq:def-Vs-norm']}, coincides with $\left( | v |_{H^1(\Omega_{1})}^2 + | v |_{V^s(\Omega_{2} | {\mathbb R}^n)}^2 \right)^{1/2} .$ We repeat the same three steps as in §\ref{['subsec_Euler-Lagrange_I']}. Step 1. Taking $v\in{\mathcal{D}}(\Omega_{1})$, the integral over $\Omega_{1}$ in \ref{['eq:weak-EIIv3']} is handled as before by a standard integration by parts. On the other hand, there is a contribution coming from the double integral over $Q_{\Omega_{2}}$, namely over $(\Omega_{1}\times\Omega_{2})\cup(\Omega_{2}\times\Omega_{1})$. Recalling the definition \ref{['eq:def-wNs']} of the weighted nonlocal Neumann operator, which acts as a flux density from $\Omega_{2}$ to $\Omega_{1}$, we conclude that, in $L^2(\Omega_{1})$, $-{\rm div}\, (\sl \nabla u) + \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u = f \hbox{in}\Omega_{1}.$ Step 2. Taking $v\in{\mathcal{D}}(\Omega_{2})$ in \ref{['eq:def-L']} and using \ref{['eq:identity-L-vanishes']}, one has $\frac{C(n,s)}{2}\iint_{Q_{\Omega_{2}}} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx = \int_{\Omega_{2}} f v \, dx.$ On the other hand, since $v$ vanishes on $\Omega_{2}^c$, formula \ref{['eq:def-Nsnl']} yields $\frac{C(n,s)}{2} \iint_{Q_{\Omega_{2}}} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx = \langle (-\Delta_{\sigma_{n\ell}})^s u, v \rangle.$ Thus, we obtain that in ${\mathcal{D}}'(\Omega_{2})$, and even in $L^2(\Omega_{2})$ since $f_{|\Omega_{2}}\in L^2(\Omega_{2})$, $(-\Delta_{\sigma_{n\ell}})^s u = f \hbox{in}\Omega_{2}.$ To obtain the Euler-Lagrange equations for $E_{II}$, we still need to consider $v\in{\mathcal{D}}(\Omega)$ in \ref{['eq:weak-EIIv3']} and perform the corresponding integration by parts. However, over $\Omega_{1}$, we only know that the sum $\mathcal{A} u := -{\rm div}\, (\sl \nabla u) + \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u$ belongs to $L^2(\Omega_{1})$ but not the terms separately. We shall proceed similarly to the second part of §\ref{['sec:traces-liftings']} to derive a integration by parts formula for $\mathcal{A}$ and deduce the form of its associated Neumann operator. Indeed, we let $\mathcal{J}^s(\Omega_{2}; \Omega_{1})$ be the closure of $C^\infty_c({\mathbb R}^n)$ with respect to the norm $\| u \|_{\mathcal{J}^s(\Omega_{2};\Omega_{1})} := \left( \| u \|_{V^s(\Omega_{2} | {\mathbb R}^n)}^2 + \| (-\Delta_{\sigma_{n\ell}})^s u \|_{L^2(\Omega_{2})}^2 + \| u \|_{H^1(\Omega_{1})}^2 + \| \mathcal{A} u \|_{L^2(\Omega_{1})}^2 \right)^{\frac{1}{2}}.$ We remark that, above, similarly to §\ref{['sec:domain-Nsnl']}, instead of $\| (-\Delta_{\sigma_{n\ell}})^s u \|_{L^2(\Omega_{2})}$, for the sake of the following argument we could have equivalently used $\| (-\Delta_{\sigma_{n\ell}})^s u \|_{H^{-\tau}(\Omega_{2})}$ for any $\tau \in (0,1/2)$. We point out that $\mathcal{J}^s(\Omega_{2};\Omega_{1}) \subset \mathcal{H}^s(\Omega_{2})$ (see definition \ref{['eq:def-calHt']}), and therefore for all $u \in \mathcal{J}^s(\Omega_{2};\Omega_{1})$ we have $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u \in \left( \mathcal{T}^{s}(\Omega_{2}^c)\right)'$ by \ref{['eq:def-Nsnl']}. We let $H^1_{0,{\Gamma_1}}(\Omega_{1})$ be the completion of $\{ v \in C^\infty_c(\overline{\Omega_{1}}) \colon v=0 \hbox{in a neighborhood of} {\Gamma_1}\}$ in $H^1(\Omega_{1})$. Since ${\Gamma_1}$ is a Lipschitz submanifold of $\partial\Omega_{1}$, we know that $H^1_{0,{\Gamma_1}}(\Omega_{1})$ is algebraically and topologically equal to the set of functions in $H^1(\Omega_{1})$ that have zero trace on ${\Gamma_1}$. Then, we define a bilinear operator $L \colon \mathcal{J}^s(\Omega_{2}; \Omega_{1}) \times H^1_{0,{\Gamma_1}}(\Omega_{1}) \to {\mathbb R}$, $L(u,v) := \langle \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u , v \rangle - \int_{\Omega_{1}} \mathcal{A} u v + \int_{\Omega_{1}} \sl \nabla u \cdot \nabla v .$ We claim that $L$ is well-defined. Given $u \in \mathcal{J}^s(\Omega_{2};\Omega_{1})$ and $v \in H^1_{0,{\Gamma_1}}(\Omega_{1})$, by the definitions of these spaces the only term in the formula above that we need to inspect with detail is the first one and, since $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u \in \left( \mathcal{T}^{s}(\Omega_{2}^c)\right)'$, it suffices to show that $v \in \mathcal{T}^{s}(\Omega_{2}^c)$. Because $v \in H^1_{0,{\Gamma_1}}(\Omega_{1})$, we can extend it first to $H^1_0(\Omega)$ by lifting $v|_{\Sigma}\in\widetilde{H}^{1/2}(\Sigma)$ to $\Omega_{2}$ and then by zero on $\Omega^c$. Denoting by $\tilde{E} v$ be such an extension, we find that $\tilde{E}v \in \mathcal{T}^{s}(\Omega_{2}^c)$ because $v|_{\Omega_{2}^c} \in H^s(\Omega_{2}^c)$, and $\| \tilde{E} v \|_{\mathcal{T}^{s}(\Omega_{2}^c)} \lesssim \|v \|_{H^1(\Omega_{1})}$. We remark that, according to \ref{['eq:def-Nsnl']}, the way $v$ is extended to $\Omega_{2}$ is irrelevant for the computation of the duality term $\langle \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u , v \rangle$. From the preceding discussion, it follows that the operator $L$ is bounded: $| L(u,v) |\le \| \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u \|_{\left(\mathcal{T}^{s}(\Omega_{2}^c)\right)'} \| \tilde{E} v \|_{\mathcal{T}^{s}(\Omega_{2}^c)}\quad + \| \mathcal{A} u \|_{L^2(\Omega_{1})} \| v \|_{L^2(\Omega_{1})}\quad + \| \sl \|_{L^\infty(\Omega_{1})} |u|_{H^1(\Omega_{1})} |v|_{H^1(\Omega_{1})}\lesssim \| u \|_{\mathcal{J}^s(\Omega_{2};\Omega_{1})} \|v \|_{H^1(\Omega_{1})},$ where we used the fact that $\| \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u \|_{\left(\mathcal{T}^{s}(\Omega_{2}^c)\right)'}\lesssim \left( \| u \|_{V^s(\Omega_{2} | {\mathbb R}^n)} + \| (-\Delta_{\sigma_{n\ell}})^s u \|_{L^2(\Omega_{2})} \right).$ Similarly to §\ref{['sec:domain-Nsnl']}, if $u \in C^2(\overline\Omega)$ and $v \in \mathcal{D}(\Omega_{1})$ we have, for the operator $L$ in \ref{['eq:def-L-II']}, that $L(u,v) = 0$. Indeed, we can split the operator $\mathcal{A} u = -{\rm div}\, (\sl \nabla u) + \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u$, cancel out the integrals involving $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u$, and integrate by parts the local term. Therefore, $L$ vanishes on $\mathcal{J}^s(\Omega_{2};\Omega_{1}) \times H^1_0(\Omega_{1})$. Next, we use the operator $L$ to define a local Neumann operator over $\Sigma$ for functions in $\mathcal{J}^s(\Omega_{2};\Omega_{1})$, which we shall denote by $D_\Sigma$ and which we expect to be in some space with differentiability order $-1/2$. Given $g \in \widetilde{H}^{1/2}(\Sigma)$, we lift it to $Eg \in H^1_{0,{\Gamma_1}}(\Omega_{1})$ and we can then define $D_\Sigma \colon\mathcal{J}^s(\Omega_{2};\Omega_{1}) \to (\widetilde{H}^{1/2}(\Sigma))'$ by $\langle D_\Sigma u , g \rangle := L(u,Eg).$ It is clear that this definition does not depend on the lifting operator $E$: given any two liftings of $g \in \widetilde{H}^{1/2}(\Sigma)$, $E_1 g, E_2 g \in H^1_{0,{\Gamma_1}}(\Omega_{1})$, we have $E_1 g - E_2 g \in H^1_0(\Omega_{1})$ and thus $L(u,E_1g - E_2 g) = 0$ for all $u \in \mathcal{J}^s(\Omega_{2};\Omega_{1})$. We are now in position to proceed with the derivation of the Euler-Lagrange equations associated to $E_{II}$. Step 3. From Step 2, we have $(-\Delta_{\sigma_{n\ell}})^s u \in L^2(\Omega_{1})$ and therefore $u \in \mathcal{J}^s(\Omega_{2};\Omega_{1})$. Thus, we take $v \in \mathcal{D}(\Omega)$ in \ref{['eq:weak-EIIv3']}, use the integration by parts formula from Lemma \ref{['lemma:ibp-snl']} and the definition \ref{['eq:def-L-II']} of the operator $L$: $\int_\Omega f v \, dx= \frac{C(n,s)}{2} \iint_{Q_{\Omega_{2}}} {\sigma_{n\ell}}(x,y) \frac{(u(x)-u(y))(v(x)-v(y))}{|x-y|^{n+2s}} dy dx + \int_{\Omega_{1}} \sl \nabla u \cdot \nabla v\, dx= \int_{\Omega_{2}} v (-\Delta_{\sigma_{n\ell}})^s u + \int_{\Omega_{1}} v \mathcal{A} u + L(u, v) .$ Because $(-\Delta_{\sigma_{n\ell}})^s u = f$ in $\Omega_{2}$ and $\mathcal{A} u = f$ in $\Omega_{1}$, we get $L(u, v) = 0$, and we recover the interface condition $\langle D_\Sigma u , v \rangle_{ \widetilde{H}^{1/2}(\Sigma)} = 0 \quad \forall v \in \mathcal{D}(\Omega).$ To conclude, the solution $u$ to our problem with the second energy is governed by: find $u\in \mathcal{H}_{II}$ such that $\left\lbrace -{\rm div}\, (\sl \nabla u) + \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u = f\hbox{in} \Omega_{1},(-\Delta_{\sigma_{n\ell}})^s u = f\hbox{in} \Omega_{2},D_\Sigma u = 0\hbox{on} \Sigma. \right.$ Let $u \in \mathcal{J}^s(\Omega_{2};\Omega_{1})$. Since $\mathcal{A} u \in L^2(\Omega_{1})$, if either ${\rm div}\, (\sl \nabla u)$ or $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u$ belong to $L^2(\Omega_{1})$, then the other one belongs to $L^2(\Omega_{1})$ as well. In such a case, $D_\Sigma u$ coincides with the (outer) normal derivative of $u$ on $\Sigma$ from $\Omega_{1}$. Indeed, splitting $\mathcal{A} u$, canceling out the integrals involving $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u$, and integrating by parts the local term we reach $\langle D_\Sigma u , g \rangle_{ \widetilde{H}^{1/2}(\Sigma)} = L(u,Eg) = \langle (\sl \frac{\partial u}{\partial {\mathbf n}})_{|\Omega_{1}}, E g \rangle_{H^{1/2}(\partial\Omega_{1})} = \langle (\sl \frac{\partial u}{\partial {\mathbf n}})_{|\Omega_{1}}, g \rangle_{ \widetilde{H}^{1/2}(\Sigma)}.$ We therefore have for all $u \in \mathcal{J}^s(\Omega_{2};\Omega_{1})$ $D_\Sigma u = (\sl \frac{\partial u}{\partial {\mathbf n}})_{|\Omega_{1}}\hbox{provided that} {\rm div}\, (\sl \nabla u)\in L^2(\Omega_{1}) \hbox{or} \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u\in L^2(\Omega_{1}).$ If $s< 1/2$ and $f|_{\Omega_{2}}$ is smoother than just $L^2$, then we actually expect the solutions to exhibit a standard homogeneous Neumann condition on the interface. Indeed, in that case one can show that $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u |_{\Omega_{1}} \in L^2(\Omega_{1})$. This, in turn, allows us to apply known results for local operators with mixed boundary conditions to prove the regularity of solutions. We postpone a full discussion of this topic to Section \ref{['sec:regularity-small-s']}. In general, we must point out that we expect neither ${\rm div}\, (\sl \nabla u)$ nor $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u$ to belong to $L^2(\Omega_{1})$ and, therefore, we do not expect the boundary condition $D_\Sigma u = 0$ on $\Sigma$ to correspond to a classical homogeneous Neumann condition. This is consistent with the fact that, as $s \to 1$, the system \ref{['eq:coupled_problemIII_NEW']} should formally converge to \ref{['eq:local-transmission']}. In this section, we shall focus on the energy $E_{II}$. We address the regularity of the minimizers to \ref{['eq:energy-III']} over $\mathcal{H}_{II}$. We shall make use of the Euler-Lagrange equations \ref{['eq:coupled_problemIII_NEW']} and some known results on regularity of solutions to local and nonlocal problems on Lipschitz domains in an iterative fashion. In first place, we review these results, and afterwards combine them to deduce the regularity of solutions under some additional assumptions on the domain configuration. Since \ref{['eq:coupled_problemIII_NEW']} involves the interface operator $D_\Sigma$, which may not coincide with a classical normal derivative, we first consider the case in which $\Sigma = \emptyset$ and therefore one recovers two isolated Dirichlet problems. Afterwards, we show that if $s<1/2$ and $f$ is sufficiently smooth we can actually prove that $D_\Sigma$ corresponds to a classical Neumann operator and derive regularity of solutions in that case. We present the results in this section for specific settings of particular interest to us. However, several variants are possible with minor adjustments to certain parameters, such as the type of domain (polygonal or arbitrary Lipschitz) or the regularity of the data. This is especially true for our choice of Proposition \ref{['prop:local-regularity-mixed']}, which could be replaced with any suitable regularity estimate applicable to problems with mixed Dirichlet/Neumann boundary conditions. Here, we briefly discuss some global regularity estimates in the Besov scale for solutions of PDE on Lipschitz domains that we shall need in the sequel. For the sake of completeness, we collected some preliminary material on Besov spaces in Appendix \ref{['sec:Besov']}. There is a vast literature concerning regularity of local, linear, elliptic PDE on bounded, non-smooth or piecewise smooth domains, cf. for example JeKe95Savare98Eb99EbFr99Jo99DiKaRe15ScSz22 and the books GrisvardCostabel_book. In this work, depending on the configuration of the domain, we need to resort to either estimates with homogeneous Dirichlet conditions (in the case $\Sigma = \emptyset$) or with mixed Dirichlet and Neumann conditions with regions that may or may not intersect, depending on whether $\overline\Sigma \cap \overline{\partial \Omega}$ is empty or not. In particular, the way in which regions with Neumann and Dirichlet boundary conditions intersect critically affects the resulting regularity of solutions KaMa07. Besides the condition $\sigma_{\max}\ge\sl \ge\sigma_{\min}>0$ that we introduced in Section \ref{['sec:l-nl-coupling']}, we assume the diffusivity $\sl$ to be Lipschitz continuous, $|\sl(x) - \sl(y) | \lesssim |x-y| , \quad \forall x, y \in \Omega_{1}.$ For problems with homogeneous Dirichlet boundary conditions, we have the following result Savare98 (see also JeKe95 for related results in spaces with different integrability indices). Since for the sake of approximation we can only exploit minimal regularity, we shall conform ourselves with results of such a type. Let $\Lambda$ be a bounded Lipschitz domain, $\sl$ satisfy \ref{['eq:Lipschitz-sigma']}, let $f \in H^{-1}(\Lambda)$ be given and $u \in H^1_0(\Lambda)$ be the weak solution to the homogeneous Dirichlet problem $\left\lbrace -{\rm div}\, (\sl \nabla u) = f\hbox{in} \Lambda,u = 0\hbox{on} \partial\Lambda. \right.$ If $f\in B^{-1/2}_{2,1}(\Lambda)$, then $u \in \dot B^{3/2}_{2,\infty}(\Lambda)$, with $\| u \|_{\dot{B}^{3/2}_{2,\infty}(\Lambda)} \lesssim \| f \|_{B^{-1/2}_{2,1}(\Lambda)}.$ If $f\in H^{-1+\frac{\theta}{2}}(\Lambda)$ for some $\theta \in (0,1)$, then $u \in \widetilde{H}^{1+ \frac{\theta}{2}}(\Lambda)$ with $\| u \|_{\widetilde{H}^{1+ \frac{\theta}{2}}(\Lambda)} \lesssim \| f \|_{H^{-1+ \frac{\theta}{2}}(\Lambda)}.$ We next list some results regarding problems with mixed boundary conditions. Since our goal is to apply the resulting regularity estimates to derive a priori approximation rates for finite element discretizations of problem \ref{['eq:coupled_problemIII_NEW']} on two-dimensional polygonal domains, we will only state results in unweighted, $L^2$-Sobolev spaces and for problems with a constant diffusivity, cf. Gr76. Other related and important results in that could also be applied to our problem include regularity for curved polyhedra in 3d Da92, or classical estimates on weighted spaces Ko67. Let $\Lambda \subset {\mathbb R}^2$ be a bounded polygonal domain whose boundary is the union of open sets $\Gamma_D$, $\Gamma_N$, with $\Gamma_D \cap \Gamma_N = \emptyset$, $\Gamma_D \neq \emptyset$, and $\overline \Gamma_D \cap \overline\Gamma_N$ consisting of a finite set of points. Let $f \colon \Lambda \to {\mathbb R}$ be given and $u \in H^1(\Lambda)$ be the weak solution to the homogeneous Dirichlet problem $\left\lbrace - \Delta u = f\hbox{in} \Lambda,u = 0\hbox{on} \Gamma_D,\frac{\partial u}{\partial {\mathbf n}} = 0\hbox{on} \Gamma_N. \right.$ Let $\mathcal{C} \subset \partial \Omega$ be the (finite) set of intersecting points between $\overline\Gamma_D$ and $\overline\Gamma_N$ and corner points in $\Gamma_D$ and $\Gamma_N$. For all $c \in \mathcal{C}$, let $\omega_c \in (0, 2 \pi)$ be the opening of the ($\Lambda$-interior) angle at $c$ and let $\gamma_c = \left\lbrace 1 + \frac{\pi}{2\omega_c}\hbox{if} c \in \overline \Gamma_D \cap \overline\Gamma_N ,1+ \frac{\pi}{\omega_c}\hbox{if} c \in \Gamma_D \cup \Gamma_N . \right.$ Then, for $\gamma = \gamma(\Lambda) = \min\{\gamma_c : \, c \in \mathcal{C} \}>1,$ if $f \in H^{\gamma-2}(\Lambda)$ we have $u \in H^\gamma(\Lambda)$ and $\| u \|_{H^\gamma(\Lambda)} \lesssim \| f \|_{H^{\gamma-2}(\Lambda)}.$ Regularity results concerning nonlocal operators such as $(-\Delta_{\sigma_{n\ell}})^s$ are much scarcer than their local counterparts. Unlike what is needed for the local operator, independently of whether $\Sigma = \emptyset$ or not, for the purposes of our paper we only require results involving Dirichlet boundary conditions, although of non-homogeneous type. At this point, we need to make some assumption on the diffusivity ${\sigma_{n\ell}}$. Following Savare98 and BoLiNo23, to obtain extra regularity of solutions we require some assumptions about the smoothness of the diffusivities. In the setting of BoLiNo23, we have $G(x,y,\rho) = \frac{C(n,s)}{4} {\sigma_{n\ell}}(x,y) \rho^2$. Therefore, to apply the results in that paper, we only need to assume some Hölder continuity on ${\sigma_{n\ell}}$. More precisely, we assume that there exists $\beta \in (0,1]$ such that $| {\sigma_{n\ell}}(x,y) - {\sigma_{n\ell}}(x',y') | \lesssim \left(|x-x'|^\beta + |y-y'|^\beta \right).$ With that, we can apply BoLiNo23 to deduce regularity estimates for problems with homogeneous Dirichlet conditions. We briefly state a result in this regard. We refer to BoNo23 and BoLiNo23. Let $s \in (0,1)$. Assume ${\sigma_{n\ell}}$ satisfies \ref{['eq:Holder-sigma']}, and let $f\in H^{-s}(\Lambda)$ be given. Let $u \in \widetilde{H}^s(\Lambda)$ be the weak solution to the homogeneous Dirichlet problem $\left\lbrace (-\Delta_{\sigma_{n\ell}})^s u = f\hbox{in} \Lambda,u = 0\hbox{in} \Lambda^c. \right.$ If $f\in B^{-s+\beta/2}_{2,1}(\Lambda)$, then $u \in \dot B^{s+\beta/2}_{2,\infty}(\Lambda)$, with $\| u \|_{\dot{B}^{s+ \frac{\beta}{2}}_{2,\infty}(\Lambda)} \lesssim \| f \|_{B^{-s+\frac{\beta}{2}}_{2,1}(\Lambda)}.$ If $f\in H^{-s+\lambda \frac{\beta}{2}}(\Lambda)$ for some $\lambda \in (0,1)$, then $u \in \widetilde{H}^{s+ \lambda \frac{\beta}{2}}(\Lambda)$ with $\| u \|_{\widetilde{H}^{s+ \lambda \frac{\beta}{2}}(\Lambda)} \lesssim \| f \|_{H^{-s+\lambda \frac{\beta}{2}}(\Lambda)}.$ The regularity estimate \ref{['eq:bes-regularity-max']} is the maximal one can expect even for a smooth right hand side $f$ and a smooth domain. As a matter of fact, the reduced regularity of solutions is caused by a rough boundary behavior, as illustrated by a simple example: $\Lambda = B(0,1)\subset{\mathbb R}^n, \ f \equiv 1, \ {\sigma_{n\ell}} \equiv 1 \ \Rightarrow \ u(x) = c(n,s) (1-|x|^2)^s_+,$ where $c(n,s)= 2^{-2s}\Gamma(\frac{n}{2})\left(\Gamma(\frac{n+2s}{2})\Gamma(1+s)\right)^{-1}$. Regularity of solutions to problem \ref{['eq:Dirichlet-FL']} with a constant diffusivity has been thoroughly studied in recent years; in particular, reference RosOtonSerra obtains sharp estimates on the boundary behavior of solutions, that leads to weighted Hölder regularity estimates. As a conclusion, the boundary expansion $u(x) \simeq d(x, \partial\Lambda)^s \varphi(x),$ where $\varphi$ is smooth, turns out to be generic. To prove regularity results of solutions of problems with non-homogeneous exterior conditions, we only need to exploit the mapping properties of the operator $(-\Delta_{\sigma_{n\ell}})^s$. For completeness, we provide a proof of the following lemma in detail. Let $s \in (0,1/2)$, $\alpha \in (0, 1-2s)$, and the diffusivity ${\sigma_{n\ell}}$ satisfy \ref{['eq:Holder-sigma']} with $\beta \ge \alpha$. Then, the weighted fractional Laplacian \ref{['eq:def-wifl']} is a bounded operator from $H^{\alpha+2s} ({\mathbb R}^n)$ to $B^{\alpha}_{2,\infty}({\mathbb R}^n).$ See Appendix \ref{['sec:mapping-Laplacian']}. We state the following theorem in a fashion that we will use in Section \ref{['sec:regularity-small-s']}. Let $\Lambda$ be a bounded Lipschitz domain, $s \in (0,1/2)$, $g \in H^{\alpha+2s}(\Lambda^c)$ for some $\alpha \in (0, 1-2s)$, ${\sigma_{n\ell}}$ satisfy \ref{['eq:Holder-sigma']} with $\beta \ge \alpha$, $f\colon \Lambda \to {\mathbb R}$, and $u \in H^s({\mathbb R}^n)$ be the weak solution to the non-homogeneous Dirichlet problem $\left\lbrace (-\Delta_{\sigma_{n\ell}})^s u = f\hbox{in} \Lambda,u = g\hbox{in} \Lambda^c. \right.$ If $f\in H^{-s+\theta \frac{\beta}{2}}(\Lambda)$ for some $\theta \in (0,1)$, then $u \in H^{s+\min\{ \alpha+s-\epsilon, \theta\frac{\beta}{2} \}}({\mathbb R}^n)$ for all $\epsilon > 0$ and $\| u \|_{H^{s+\min\{ \alpha+s-\epsilon, \theta\frac{\beta}{2} \}}({\mathbb R}^n)} \lesssim \| f \|_{H^{-s+\theta\frac{\beta}{2}}(\Lambda)} + \| g \|_{H^{\alpha+2s} (\Lambda^c)}.$ The hidden constant above behaves as $\epsilon^{-1/2}$ if $\alpha+s \le \theta\frac{\beta}{2}$ and as $(\alpha+s-\theta\frac{\beta}{2})^{-1}$ if $\alpha +s > \theta \frac{\beta}{2}$. As $\Lambda$ is a bounded Lipschitz domain and $g \in H^{\alpha+2s}(\Lambda^c)$, we extend it to $G \in H^{\alpha+2s}({\mathbb R}^n)$. By using Lemma \ref{['lemma:mapping-Laplacian']}, we have $(-\Delta_{\sigma_{n\ell}})^s G \in B^{\alpha}_{2,\infty}({\mathbb R}^n)$ and then the function $\tilde{u} = u - G$ is a weak solution of the homogeneous Dirichlet problem $\left\lbrace (-\Delta_{\sigma_{n\ell}})^s \tilde{u} = \tilde{f}\hbox{in} \Lambda,\tilde{u} = 0\hbox{in} \Lambda^c, \right.$ with $\tilde{f} = f - (-\Delta_{\sigma_{n\ell}})^s G \in B^{\alpha}_{2,\infty}(\Lambda)$. In case $\alpha \le -s + \theta \frac{\beta}{2}$, we use the embedding $B^{\alpha}_{2,\infty}(\Lambda) \subset H^{\alpha -\epsilon}(\Lambda)$ for some arbitrary $\epsilon > 0$ (see \ref{['eq:Besov-Sobolev-emb']}), that gives rise to $\| \tilde{f} \|_{H^{\alpha -\epsilon}(\Lambda)} \lesssim \| f \|_{H^{-s + \theta \frac{\beta}{2}}(\Lambda)} + \epsilon^{-1/2} \| (-\Delta_{\sigma_{n\ell}})^s G \|_{B^{\alpha}_{2,\infty}(\Lambda)},$ and use \ref{['eq:sob-regularity-intermediate']} with $\lambda = \frac{2(\alpha+s-\epsilon)}{\beta} \in (0,1)$ to obtain $\| \tilde{u} \|_{\widetilde{H}^{\alpha +2s - \epsilon}(\Lambda)} \lesssim \| f \|_{H^{-s + \theta \frac{\beta}{2}}(\Lambda)} + \epsilon^{-1/2} \| g \|_{H^{\alpha+2s} (\Lambda^c)} .$ For $\alpha > -s + \theta \frac{\beta}{2}$ we again resort to \ref{['eq:Besov-Sobolev-emb']}, exploit the embedding $B^{\alpha}_{2,\infty}(\Lambda) \subset H^{-s+\theta\frac{\beta}{2}}(\Lambda)$ with a constant $\mathcal{O}(\alpha+s-\theta\frac{\beta}{2})$, and use \ref{['eq:sob-regularity-intermediate']} with $\lambda = \theta$: $\| \tilde{u} \|_{\widetilde{H}^{s+\theta \frac{\beta}{2}}(\Lambda)} \lesssim \| \tilde{f} \|_{H^{-s+\theta\frac{\beta}{2}}(\Lambda)} \lesssim \| f \|_{H^{-s+\theta\frac{\beta}{2}}(\Lambda)} + \frac{ \| g \|_{H^{\alpha+2s} (\Lambda^c)}}{\alpha+s-\theta\frac{\beta}{2}} .$ We proved Lemma \ref{['lemma:mapping-Laplacian']} and Theorem \ref{['thm:reg-nl-s-small']} for $s$ in the range $(0,1/2)$. The main reason behind this limitation becomes apparent upon inspection of our argument in Appendix \ref{['sec:mapping-Laplacian']}: we used first differences. For sufficiently smooth diffusivities, we expect the weighted fractional Laplacian to be an operator of order $2s$, but if $2s > 1$ then one cannot derive such a regularity by using first differences. If $\alpha \in (0, 2-2s)$ and ${\sigma_{n\ell}}$ is locally $C^\beta$ near the diagonal $\{x = y\}$ for some $\beta \ge \alpha$, then we expect $(-\Delta_{\sigma_{n\ell}})^s \colon H^{\alpha+2s} ({\mathbb R}^n) \to B^{\alpha}_{2,\infty}({\mathbb R}^n)$ to be a bounded operator. However, we emphasize that we are interested in applying Theorem \ref{['thm:reg-nl-s-small']} in the case where $\Sigma \neq \emptyset$ and the operator $D_\Sigma$ coincides with the classical Neumann operator; as we already anticipated, and show in Section \ref{['sec:regularity-small-s']} below, we expect this only to happen provided $s<1/2$. From our discussion in Section \ref{['subsec_Euler-Lagrange_III']}, and in particular from the resulting Euler-Lagrange equations \ref{['eq:coupled_problemIII_NEW']}, we note that the problem simplifies notably in the case $\overline\Omega_{1} \cap \overline\Omega_{2} = \emptyset$. Indeed, in such a case, for all $x \in \Omega_{1}$ and $y \in \Omega_{2}$ we have $|x-y| \ge d(\Omega_{1}, \Omega_{2}) > 0$ and thus we have the pointwise bound $|\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u (x)|\le \frac{C(n,s) \, \| {\sigma_{n\ell}} \|_{L^\infty(\Omega_{1} \times \Omega_{2})}}{d(\Omega_{1}, \Omega_{2})^{n+2s}} \int_{\Omega_{2}} |u(x)-u(y)| dy\lesssim \left[ |u(x)| + \| u \|_{L^2(\Omega_{2})} \right]$ for all $x \in \Omega_{1}$, with a hidden constant depending on $n,s, {\sigma_{n\ell}}, d(\Omega_{1}, \Omega_{2}),$ and $|\Omega_{2}|$. Therefore, because $u \in L^2(\Omega)$ with $\| u \|_{L^2(\Omega)}\lesssim \| f \|_{L^2(\Omega)}$ --this is because the energy space satisfies $\mathcal{H}_{II}\subset L^2(\Omega)$--, we immediately deduce $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u \in L^2(\Omega_{1}) , \quad \hbox{with} \quad \| \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u \|_{L^2(\Omega_{1})} \lesssim \| f \|_{L^2(\Omega)}.$ Therefore, in this setting, in order to prove regularity of solutions within $\Omega_{1}$ it suffices to apply Proposition \ref{['prop:local-regularity-dir']}. Additionally, in the subdomain $\Omega_{2}$ we can exploit the fact that $u|_{\Omega_{2}^c}$ vanishes on a neighborhood of $\Omega_{2}$. Namely, if we denote by $U$ the zero-extension of $u|_{\Omega_{2}^c}$ to $\Omega_{2}$, for $x \in \Omega_{2}$ we have $(-\Delta_{\sigma_{n\ell}})^s U (x)= C(n,s) \int_{{\mathbb R}^n} {\sigma_{n\ell}}(x,y) \frac{U(x) - U(y)}{|x-y|^{n+2s}} \, dy= - C(n,s) \int_{\Omega_{1}} {\sigma_{n\ell}}(x,y) \frac{U(y)}{|x-y|^{n+2s}} \, dy$ and, because $d(\Omega_{1}, \Omega_{2}) > 0$, if ${\sigma_{n\ell}}$ satisfies \ref{['eq:Holder-sigma']} we deduce $(-\Delta_{\sigma_{n\ell}})^s U |_{\Omega_{2}} \in C^\beta(\overline\Omega_{2}) \subset H^{\beta-\epsilon}(\Omega_{2})$ for all $\epsilon > 0$. We can therefore mimic the proof of Theorem \ref{['thm:reg-nl-s-small']} with this modification and arbitrary $s \in (0,1)$, and obtain enhanced regularity estimates. We summarize our discussion in the following result. Let $s \in (0,1)$. Assume $\overline\Omega_{1} \cap \overline \Omega_{2} = \emptyset$, $\beta \in [0,1)$ with $\beta \ge 2s - 1$, $f \in B^{-s+\beta/2}_{2,1}(\Omega)$, and let $u \in \mathcal{H}_{II}$ be the minimizer of \ref{['eq:energy-III']}. Then, $u|_{\Omega_{1}} \in \dot B^{3/2}_{2,\infty}(\Omega_{1})$ and $u|_{\Omega_{2}} \in \dot B^{s+\beta/2}_{2,\infty}(\Omega_{2})$. Let $u \in \mathcal{H}_{II}$ be the minimizer of \ref{['eq:energy-III']}. According to our discussion above, on $\Omega_{1}$ we have $-{\rm div}\, (\sl \nabla u) = f - \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u \in B^{-1/2}_{2,1}(\Omega_{1}), \quad u |_{\partial \Omega_{1}} = 0,$ and applying Proposition \ref{['prop:local-regularity-dir']} we deduce $u|_{\Omega_{1}} \in \dot B^{3/2}_{2,\infty}(\Omega_{1})$. To prove the claimed regularity of $u|_{\Omega_{2}}$, we point out that the function $\tilde{u} = u - U$ satisfies the homogeneous Dirichlet problem $(-\Delta_{\sigma_{n\ell}})^s \tilde{u} = f - (-\Delta_{\sigma_{n\ell}})^s U \in B^{-s+\beta/2}_{2,1}(\Omega_{2}), \quad \tilde{u} |_{\Omega_{2}^c} = 0,$ and we can apply \ref{['eq:bes-regularity-max']}. The assumption $\beta \ge 2s - 1$ in the previous proposition is by no means fundamental. We included it for the clarity of the statement. By the same argument, if $s>1/2$ one can also prove regularity in case $\beta < 2s - 1$, only that in such a case one would have a reduced regularity over $\Omega_{1}$. Similarly, we can also apply interpolation arguments to deduce regularity estimates in case $f$ is less regular than $B^{-s+\beta/2}_{2,1}(\Omega)$. We now assume that $s \in (0,1/2)$, $\alpha \in (2s, s+1/2)$, and $f \in H^{\alpha - 2s}(\Omega_{2}) \subset L^2(\Omega_{2})$. As anticipated in Remark \ref{['rem:homogeneous-Neumann']}, because $s<1/2$ we can actually show that the solution $u$ satisfies a homogeneous Neumann boundary condition on the interface from $\Omega_{1}$. Let $\Omega$ be a bounded Lipschitz domain, ${\sigma_{n\ell}}$ satisfy \ref{['eq:Holder-sigma']}, $u$ be a minimizer of the energy $E_{II}$ with $s \in (0, 1/2)$ and $f \in H^{-s+\theta\frac{\beta}{2}}(\Omega)$ for some $\theta \in (0,1)$ and satisfying $s < \theta\frac{\beta}{2}$. Then, $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u |_{\Omega_{1}} \in L^2(\Omega_{1})$. The well-posedness of problem \ref{['eq:coupled_problemIII_NEW']} implies $u|_{\Omega_{2}^c} \in H^1(\Omega_{2}^c)$, and because $u$ vanishes on $\Omega^c$, we have $u|_{\Omega_{2}^c} \in H^{1-\epsilon}(\Omega_{2}^c)$ for all $\epsilon > 0$. We can therefore apply Theorem \ref{['thm:reg-nl-s-small']} with $\alpha = 1 - 2s - \epsilon$; because $s < 1/2$ we can assume $\epsilon$ is small enough so that $\min\{ 1- s - \epsilon, \theta\frac{\beta}{2} \} = \theta\frac{\beta}{2}$ and therefore we obtain the global regularity $u \in H^{s+\theta\frac{\beta}{2}}({\mathbb R}^n)$. Let us define $\omega \colon \Omega_{1} \to {\mathbb R}$, $\omega (x) = \left(\int_{\Omega_{2}} \frac{1}{|x-y|^{n+2s- \theta\beta}} \, dy \right)^{1/2}.$ Then, because $\theta\frac{\beta}{2} > s$ and $\Omega_{1}$, $\Omega_{2}$ are bounded, an integration in polar coordinates allows us to deduce that $\omega \in L^\infty(\Omega_{1})$ (with $\| \omega \|_{L^\infty(\Omega_{1})} \le C( n, s, \alpha, d(\Omega_{1}), d(\Omega_{2}) )$). We can therefore compute, for all $x \in \Omega_{1}$ and by applying Hölder's inequality, $|\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u (x)|\le C(n,s) \, \| {\sigma_{n\ell}} \|_{L^\infty(\Omega_{1} \times \Omega_{2})} \int_{\Omega_{2}} \frac{|u(x)-u(y)|}{|x-y|^{n+2s}} \, dy\le C(n,s) \, \| {\sigma_{n\ell}} \|_{L^\infty(\Omega_{1} \times \Omega_{2})} \, \| \omega \|_{L^\infty(\Omega_{1})} \left(\int_{\Omega_{2}} \frac{|u(x)-u(y)|^2}{|x-y|^{n+2s+\theta\beta}} \, dy \right)^{1/2}.$ Squaring both sides above and integrating on $\Omega_{1}$ yields $\| \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u \|_{L^2(\Omega_{1})} \lesssim \left(\int_{\Omega_{1}} \int_{\Omega_{2}} \frac{|u(x)-u(y)|^2}{|x-y|^{n+2s+\theta\beta}} \, dy \, dx \right)^{1/2}.$ The last integral above can be bounded by $\| u \|_{H^{s+\theta\frac{\beta}{2}}({\mathbb R}^n)}$ and we conclude that $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u \in L^2(\Omega_{1})$. Under the assumptions of the previous lemma, by (\ref{['eq:normal-derivative']}) we deduce that, on $\Omega_{1}$, $u$ solves the problem $\left\lbrace -{\rm div}\, (\sl \nabla u) = f - \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u\hbox{in} \Omega_{1},(\sl \frac{\partial u}{\partial {\mathbf n}})|_{\Omega_{1}} = 0\hbox{on} \Sigma,u = 0\hbox{on} {\Gamma_1} = \partial\Omega_{1} \setminus \Sigma. \right.$ We have $f - \mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}} u \in L^2(\Omega_{1})$. Next, we simplify our discussion by restricting $\Omega_{1}$ to be a two-dimensional, polygonal domain so that Proposition \ref{['prop:local-regularity-mixed']} can be applied. Let $\Omega \subset {\mathbb R}^2$ be a bounded polygonal domain such that $\Omega_{1}$ satisfies the assumptions in Proposition \ref{['prop:local-regularity-mixed']}. Let $f \in H^{-s+\theta\frac{\beta}{2}}(\Omega)$ for some $\theta \in (0,1)$ and satisfying $s < \theta\frac{\beta}{2}$, and $u$ be a minimizer of the energy $E_{II}$ with $s<1/2$. Then, $u|_{\Omega_{1}} \in H^{\gamma}(\Omega_{1})$ for $\gamma = \min\{\gamma(\Omega_{1}),2\}$ and $u \in H^{s+\theta\frac{\beta}{2}}({\mathbb R}^2)$. Let $\gamma = \min\{\gamma(\Omega_{1}),2\} \in (1,2]$ be as in Proposition \ref{['prop:local-regularity-mixed']}. We have $f|_{\Omega_{1}} \in L^2(\Omega_{1})$ and therefore $u|_{\Omega_{1}} \in H^{\gamma}(\Omega_{1})$. This readily implies $u|_{\Omega_{2}^c} \in H^{\gamma}(\Omega_{2}^c)$ and $u|_{\Omega_{2}^c} \in H^{1-\epsilon}(\Omega_{2}^c)$ for all $\epsilon > 0$. We can apply Theorem \ref{['thm:reg-nl-s-small']} to conclude that $u \in H^{s+\theta\frac{\beta}{2}}({\mathbb R}^2)$. In this subsection, we have performed two different simplifications. On the one hand, we have restricted the analysis to problems on two spatial dimensions, and on polygonal domains. As we already stated, we made this assumption with the only goal of simplifying our analysis, and the proof we made could be easily adapted for other configurations in which different type of regularity estimates on $\Omega_{1}$ are available. On the other hand, our analysis only covered the case $s<1/2$. The reason for this assumption is twofold. First, in our derivation of the strong form of the Euler-Lagrange equations we were able to show the transmission condition $D_\Sigma = 0$, and we actually expect the operator $D_\Sigma$ to coincide with a Neumann derivative only when $s<1/2$; it should be noted that, in the limit $s \to 1$ our problem recovers a standard transmission problem and therefore the interface condition we expect is like in \ref{['eq:local-transmission']}. Second, the mapping properties of the operator $\mathcal{N}_{s,{\sigma_{n\ell}}}^{\Omega_{2}}$ are not clear when $s \ge 1/2$: while in principle we expect it to be an operator of order $2s \ge 1$, the lack of second differences in its definition difficults the analysis. In particular, we are not aware of a result like Theorem \ref{['thm:trace-extension']} for $t > 1$ and therefore we cannot give a definition like \ref{['eq:def-Nsnl']} for functions in $\mathcal{H}^t(\Omega_{2})$ with $t > 1$. Nevertheless, a direct treatment of \ref{['eq:coupled_problemIII_NEW']} as a transmission problem might allow to obtain regularity estimates in the whole range $s \in (0,1)$. In this section, we consider a finite element discretization of both proposed energies. We discuss their approximation capabilities and outline the main ingredients required in their analysis. We also include several numerical experiments, using the code developed in AcBeBo17 as a starting point, that illustrate the qualitative behavior of the proposed models as well as some of their most relevant features. We consider discretizations of the problems \ref{['eq:coupled_problemI_NEW']} and \ref{['eq:coupled_problemIII_NEW']} by means of the finite element method with piecewise linear continuous functions. Let $h_0 > 0$; for $h \in (0, h_0]$, we let $\mathcal{T}_h$ denote a mesh of $\Omega$, i.e., $\mathcal{T}_h = \{T\}$ is a partition of $\Omega$ into simplices $T$ of diameter $h_T$ and $h = \max_{T \in {\mathcal{T}_h}} h_T$. Additionally, we assume that $\mathcal{T}_h$ meshes exactly both $\Omega_{1}$ and $\Omega_{2}$, and that the family $\{{\mathcal{T}_h} \}_{h>0}$ is shape-regular, namely, $\sigma := \sup_{h>0} \max_{T \in {\mathcal{T}_h}} \frac{h_T}{\rho_T} <\infty,$ where $\rho_T$ is the diameter of the largest ball inscribed in $T$. We take simplices to be closed sets. We point out that our only assumption on meshes is shape-regularity; in particular, we do not assume any uniformity. In some cases, we choose meshes that are suitably graded towards the interface $\Sigma$ to better capture the solution behavior near that region. Let $\overline{\mathcal{N}_h}$ be the set of vertices of ${\mathcal{T}_h}$, $\mathcal{N}_h$ be the set of interior vertices, $N = \#\mathcal{N}_h$ and $\{ \varphi_i \}_{i=1}^N$ be the standard piecewise linear Lagrangian basis, with ${\varphi}_i$ associated to the node $\texttt{x}_i \in \mathcal{N}_h$. With this notation, the set of discrete functions is $\widetilde{\mathbb{V}}_h := \left\{v_h: {\mathbb R}^n \to {\mathbb R} \colon v_h \in C({\mathbb R}^n), \ v_h = \sum_{i=1}^N v_i \varphi_i \right\},$ where $v_h$ is trivially extended by zero outside $\Omega$. It is clear that $\widetilde{\mathbb{V}}_h \subset \widetilde{H}^s(\Omega)$ for all $s \in (0,1)$. Therefore, we consider conforming finite element discretizations and seek $u_h \in \widetilde{\mathbb{V}}_h$ such that either $\int_{\Omega_{1}} \sl \nabla u_h \cdot \nabla v_h\, dx + \frac{ C(n,s) }{2} \iint_{\Omega_{2}^2} {\sigma_{n\ell}}(x,y) \frac{(u_h(x)-u_h(y))(v_h(x)-v_h(y))}{|x-y|^{n+2s}} dy dx= \int_\Omega f v_h \, dx, \qquad \forall v_h \in \widetilde{\mathbb{V}}_h$ or $\int_{\Omega_{1}} \sl \nabla u_h \cdot \nabla v_h\, dx + \frac{C(n,s)}{2} \iint_{Q_{\Omega_{2}}} {\sigma_{n\ell}}(x,y) \frac{(u_h(x)-u_h(y))(v_h(x)-v_h(y))}{|x-y|^{n+2s}} dy dx= \int_\Omega f v_h \, dx, \quad \forall v_h \in \widetilde{\mathbb{V}}_h$ hold. These equations are simply discrete counterparts to \ref{['eq:Euler-Lagrange-E_I']} and \ref{['eq:weak-EIIv3']}, respectively. Clearly, $u_h$ solves \ref{['eq:Euler-Lagrange-E_I-discrete']} (resp. \ref{['eq:weak-EIIv3-discrete']}) if and only if it is the minimizer of the restriction of the convex functional $E_I$ from \ref{['eq:energy-I']} (resp. $E_{II}$ from \ref{['eq:energy-III']}) over the linear space $\widetilde{\mathbb{V}}_h$; existence of discrete solutions follows immediately. Moreover, if we let $v_h = u_h$ in either \ref{['eq:Euler-Lagrange-E_I-discrete']} or \ref{['eq:weak-EIIv3-discrete']}, then we immediately obtain discrete stability bounds given $f \in L^2(\Omega)$ (or even in suitable negative-order Sobolev spaces). Furthermore, as $\widetilde{\mathbb{V}}_h \subset \mathcal{H}_{I}$ (resp. $\widetilde{\mathbb{V}}_h \subset \mathcal{H}_{II}$) we trivially have a Galerkin orthogonality property and consequently a best approximation result holds. From this, one can combine suitable quasi-interpolation estimates with a density argument like RaviartThomas83 to prove the convergence of the approximations. In case one has at hand regularity estimates for the solutions, such as the ones described in Section \ref{['sec:regularity']}, convergence rates can be obtained. We do not delve into details, and only outline the main ingredients for the energy $E_{II}$. Using the fact that $\sl, {\sigma_{n\ell}}$ are bounded and uniformly positive, the best approximation result yields $| u - u_h|_{H^1(\Omega_{1})} + | u - u_h|_{V^s(\Omega_{2} | {\mathbb R}^n)} \lesssim \inf_{v_h \in \widetilde{\mathbb{V}_h}} | u - v_h|_{H^1(\Omega_{1})} + | u - v_h|_{V^s(\Omega_{2} | {\mathbb R}^n)}.$ Thus, it suffices to use a suitable interpolation operator to bound the sum in the right-hand side above. For clarity, we specifically consider the Scott-Zhang quasi-interpolation operator ScottZhang; analogous conclusions could be drawn from interpolation using e.g. the Clément Clement or Chen-Nochetto ChenNochetto operators, although some modifications would be needed to deal with the interface. Since in all our problems the solutions exhibit the minimal regularity $u \in H^\ell(\Omega)$ for some $\ell > 1/2$ for adequately smooth ${\sigma_{n\ell}}$ and $f$ (see Theorem \ref{['thm:reg-nl-s-small']} with $\alpha$ close to $1-2s$), we can follow the original construction of the interpolation operator from ScottZhang. For the following discussion, we consider an arbitrary bounded, polyhedral Lipschitz domain $\Lambda$ that is meshed exactly by a shape-regular family $\mathcal{T}_h = \{T\}$, and keep the notation introduced above. As we are considering piecewise linear Lagrange elements, for every node $\texttt{x}_i \in \overline{\mathcal{N}_h}$ we consider an $(n-1)$-simplex $K_i$ such that $\texttt{x}_i \in \partial K_i$, we design $I_h : \widetilde{H}^\ell(\Lambda) \to \widetilde{\mathbb{V}}_h$ by introducing the $L^2$-dual basis of the nodal basis over $K_i$ and, in particular, we let $\psi_i$ be the basis element such that $\int_{K_i} \varphi_i \psi_i = 1$. With this, we define $I_h v = \sum_{i \colon \texttt{x}_i \in \overline{\mathcal{N}_h}} \left( \int_{K_i} v(x) \psi_i(x) dx \right) \varphi_i.$ Even if the sum is taken over all the mesh nodes, for functions that vanish on $\partial \Lambda$ we have $I_h v |_{\partial \Omega} = 0$. For the analysis of the operator, we need the following notions: given $T \in {\mathcal{T}_h}$, we respectively define its first and second star by $S^1_T := \bigcup \left\{ T' \in {\mathcal{T}_h}: \, T \cap T' \neq \emptyset \right\}, \quad S^2_T := \bigcup \left\{ T' \in {\mathcal{T}_h}: \, S^1_T \cap T' \neq \emptyset \right\}.$ We also require a slight modification to deal with boundary-touching simplices: we define $\widetilde{S}_T^1 := \left\lbrace S_T^1\textrm{if } T \cap \partial\Omega = \emptyset,B_T\textrm{otherwise,}\right.$ where $B_T$ is a ball centered at the barycenter of $T$, with radius comparable to $h_T$, and such that $S_T^1\subset B_T$. We also consider $\widetilde{S}_T^2 := \bigcup \{ \widetilde{S}_{T'}^1 \colon \, T' \in {\mathcal{T}_h}, \, \hbox{int}(T') \cap S_T^1 \neq \emptyset \}.$ In ScottZhang, approximation properties of this operator in integer-order Sobolev spaces are studied. Combining the results therein with interpolation between Sobolev spaces, we deduce the following $\| v - I_h v \|_{L^2(T)} + h_T \| \nabla (v - I_h v) \|_{L^2(T)} \lesssim h_T^r |v|_{H^r(S^1_T)} \quad \forall v \in H^r(S^1_T), \ r \in [1, 2],$ with $\lesssim$ independent of $h$ and $T\in{\mathcal{T}_h}$. For less-regular functions, the approximation properties of the Scott-Zhang operator were first studied in Ciarlet. Here, we recall a related estimate derived in AcBo17BoNo23-2, that is tailored to the analysis of our problem: for all $v \in H^r(\widetilde{S}_T^2)$, $s \in (0,1)$, $r \in [s,2]$, it holds that $\left(\iint_{T \times \widetilde{S}_T^1} \frac{|(v-I_h v) (x) - (v-I_h v) (y)|^2}{|x-y|^{n+2s}} \, dy \, dx \right)^{1/2} \lesssim h_T^{r-s} |v|_{H^r(\widetilde{S}_T^2)}.$ To put together the local approximability estimates with respect to fractional-order seminorms, one needs to use a suitable localization. The first results along this line were obtained in Faermann2Faermann, where localization of the $H^s(\Lambda)$-seminorms are studied; a slight modification needs to be introduced to deal with $\widetilde{H}^s(\Lambda)$-seminorms, see BoNo23-2 for details. Concretely, we have $\|v\|_{\widetilde{H}^s(\Lambda)}^2 \le \frac{C(n,s)}{2} \sum_{T \in {\mathcal{T}_h}} \left[ \iint_{T \times \widetilde{S}_T^1}\! \frac{|v (x) - v (y)|^2}{|x-y |^{n+2s}} \; dy \; dx + \frac{ c(n,\sigma) \| v \|^2_{L^2(T)}}{s h_T^{2s}} \right],$ where $\sigma$ is the shape-regularity constant introduced in \ref{['eq:shape-regularity']}. Collecting the results mentioned above, we obtain the following global interpolation estimates. Let $\Lambda$ be a bounded, polyhedral Lipschitz domain, $s \in (0,1)$, and $I_h$ be the Scott-Zhang interpolation operator defined above on a shape regular family of meshes of $\Lambda$. Then, we have $|v - I_h v |_{H^1(\Lambda)} \lesssim \left(\sum_{T \in {\mathcal{T}_h}} h_T^{2(r-1)} |v|_{H^r(S^1_T)}^2 \right)^{1/2} \quad \forall v \in H^r(\Lambda), \ r \in [1, 2],$ and $|v - I_h v |_{\widetilde{H}^s(\Lambda)} \lesssim \left(\sum_{T \in {\mathcal{T}_h}} h_T^{2(r-s)} |v|_{H^r(\widetilde{S}^2_T)}^2 \right)^{1/2} \quad \forall v \in \widetilde{H}^r(\Lambda), \ r \in [s, 2].$ Under the assumption that ${\mathcal{T}_h}$ meshes exactly both $\Omega_{1}$ and $\Omega_{2}$, the result above immediately allows one to deduce convergence rates for our energy $E_{II}$ provided we have at hand some estimates on the regularity of solutions. Indeed, it is clear that \ref{['eq:SZ-H1']} yields convergence rates in $H^1(\Omega_{1})$ for the quasi-interpolation operator. Additionally, even in problems in which $\Sigma \neq \emptyset$, the seminorm \ref{['eq:def-Vs-norm']} can be trivially bounded in terms of the global $\widetilde{H}^s(\Omega)$ norm, since $Q_{\Omega_{2}} \subset Q_\Omega$. Thus, we can bound $|u - I_h u|_{V^s(\Omega_{2} | {\mathbb R}^n)} \le |u - I_h u|_{\widetilde{H}^s(\Omega)}$ and we can use \ref{['eq:SZ-Hs']}. In the settings of either Sections \ref{['sec:isolated-subdomains']} or \ref{['sec:regularity-small-s']}, the resulting convergence rates read as follows. Let $\{{\mathcal{T}_h} \}_{h>0}$ be a shape-regular, conforming family of meshes in $\Omega$, such that for every $h>0$ the meshes ${\mathcal{T}_h}$ cover exactly both $\Omega_{1}$ and $\Omega_{2}$, and let $h = \max_{T \in {\mathcal{T}_h}} h_T$. In the setting of Proposition \ref{['prop:reg-isolated']}, we have $| u - u_h|_{H^1(\Omega_{1})} + | u - u_h|_{V^s(\Omega_{2} | {\mathbb R}^n)} \lesssim h^{\frac{\beta}{2}} \| f \|_{B^{-s+\beta/2}_{2,1}(\Omega)}$ Additionally, under the same hypotheses as in Proposition \ref{['prop:regularity-small-s']}, we have $| u - u_h|_{H^1(\Omega_{1})} + | u - u_h|_{V^s(\Omega_{2} | {\mathbb R}^n)} \lesssim h^{\min\{\gamma, \, \theta\frac{\beta}{2}\}} \| f \|_{ H^{-s+\theta\frac{\beta}{2}}(\Omega)}.$ In this section we present the results to several computational experiments involving the coupled models we proposed in Section \ref{['sec:l-nl-coupling']}. We mainly focus on \ref{['eq:coupled_problemIII_NEW']}, but we also show some experiments related to \ref{['eq:coupled_problemI_NEW']}. In general, explicit solutions to the problems are not available, and therefore our purpose in this section is, rather than numerically corroborating the validity of Proposition \ref{['prop:conv-rates']}, to explore some salient features of this problem. Concretely, we investigate the behavior of solutions near the interface of the domain and provide evidence that the interface operator $D_\Sigma$ from \ref{['eq:coupled_problemIII_NEW']} coincides with the normal derivative for $s < 1/2$. We illustrate that, even when $\Omega$ consists of two disconnected subdomains, there is effective coupling between the models. We also present numerical evidence of the development of singularities on $\Omega_{1}$, caused by the way $\Sigma$ intersects $\partial \Omega$, even when $\Omega_{1}$ is a convex domain. Furthermore, we study the eigenvalues of the discrete problems for a fixed mesh size and show numerically that \ref{['eq:coupled_problemIII_NEW']} recovers a standard local transmission problem in the limit as $s \to 1$. Finally, we compare the behavior of the two models \ref{['eq:coupled_problemI_NEW']} and \ref{['eq:coupled_problemIII_NEW']} in this limit. As discussed in Section \ref{['sec:regularity-small-s']}, a relevant feature of solutions to \ref{['eq:coupled_problemIII_NEW']} is that, for $s<1/2$, there is a homogeneous Neumann interface condition from the local subdomain. In contrast, for $s \to 1$ we formally expect to recover a standard transmission condition. We illustrate this behavior with the following experiment. Let $\Omega = B(0,2)$, with $\Omega_{2} = B(0,1)$, $\Omega_{1} = B(0,2) \setminus \overline{B(0,1)}$, $\sl \equiv 1$, ${\sigma_{n\ell}} \equiv 1$, $f \equiv 1$. Formally, for $s=1$ we recover the Dirichlet problem $u \in H^1_0(\Omega)$: $-\Delta u = 1$, whose solution is $u(x_1,x_2) = \frac{(4-(x_1^2+x_2^2))_+}{4}$ and satisfies ${\frac{\partial u}{\partial {\mathbf n}}}|_{\Omega_{1}} \equiv -\frac{1}{2}$ at the interface $\Sigma = \partial B(0,1)$. We compute the normal derivative on the interface of the solutions to \ref{['eq:weak-EIIv3-discrete']} for $s = \frac{j}{100}$, $j = 1, \ldots, 99$. To better capture the solution behavior, we consider approximations on adapted meshes in the spirit of Grisvard Grisvard grading according to the distance to the interface. More precisely, we fix $h>0$ and $\mu \ge 1$ and build meshes such that the triangles have size $h_T \simeq h \,d(T,\Sigma)^{(\mu-1)/\mu}\hbox{if} T\cap \Sigma = \emptyset,h^{\mu}\hbox{if} T\cap \Sigma \neq \emptyset.$ Figure \ref{['fig:interface-two-balls']} exhibits different slices of the computed solutions, as well as the computed value of the normal derivative at the point $(1,0)$, obtained on a mesh graded as above with $\mu = 2.5$ and with $26077$ elements on $\Omega_{1}$ and $40306$ on $\Omega_{2}$. From the right panel, it is apparent that $\left|\frac{\partial}{\partial {\mathbf n}}(u_h|_{\Omega_1})(1,0)\right|$ increases with respect to $s$. While it is clear that none of the computed values equals zero, we point out that for $s=0.01$ and $s=0.5$ we obtained $\frac{\partial}{\partial {\mathbf n}}(u_h|_{\Omega_1})(1,0) \approx -6.7\times 10^{-4}$ and $\frac{\partial}{\partial {\mathbf n}}(u_h|_{\Omega_1})(1,0) \approx -9.4\times 10^{-3}$, respectively. In contrast, for $s=0.99$ we computed $\frac{\partial}{\partial {\mathbf n}}(u_h|_{\Omega_1})(1,0) \approx -0.430$; to better illustrate the convergence for $s \to 1$, we also run the experiment with $s = 0.999$ and obtained $\frac{\partial}{\partial {\mathbf n}}(u_h|_{\Omega_1})(1,0) \approx -0.494$, in very good agreement with the expected limit value $-1/2$. Finally, we point out that, in our experiments we observed that, for all $s$, $\frac{\partial}{\partial {\mathbf n}}(u_h|_{\Omega_1})(1,0)$ was monotonically increasing as the meshes were increasingly refined. Finite element solutions to the problem described in Section \ref{['sec:two-balls-experiment']}. Left panel: slices $\{ x_1 > 0 ,x_2 = 0\}$ for $s=0.1, 0.5, 0.9$. Right panel: computed values of $\frac{\partial}{\partial {\mathbf n}}(u_h|_{\Omega_1})(1,0)$ for different values of $s$. Our next set of experiments aims to illustrate the coupling between the local and nonlocal models for the energy \ref{['eq:energy-III']} even when $\Omega$ is composed of two isolated subdomains, namely $\Sigma = \emptyset$. We consider $\Omega_{1} = (-\frac{3}{4}, -\frac{1}{4}) \times (-\frac{1}{2},\frac{1}{2})$, $\Omega_{2} = (\frac{1}{4}, -\frac{3}{4}) \times (-\frac{1}{2},\frac{1}{2})$; $s = 0.5$; $\sl \equiv 1$, ${\sigma_{n\ell}} \equiv 1$; and either $f = \chi_{\Omega_{1}}$ or $f = \chi_{\Omega_{2}}$. Figure \ref{['fig:two_subdomains']} shows that finite element solutions we obtained in this setting and indicates that, even if $f$ is identically zero on one subdomain in each case, the solution does not vanish. On each panel, we use different color maps for each subdomain so that the behavior of the solution becomes clearer. Minimizers of \ref{['eq:weak-EIIv3-discrete']} in the setting of Section \ref{['sec:isolated-subdomains']}. Left panel: $f = \chi_{\Omega_{1}}$; right panel: $f = \chi_{\Omega_{2}}$. Our next experiment explores the development of algebraic singularities over $\Omega_{1}$ for minimizers to \ref{['eq:energy-III']} and their dependence on $s$. We consider $\Omega = (-1/2,1/2)^2 \setminus [-1/2,0]^2$, with $\Omega_{1} = \Omega \cap \{ x_2 > 0 \}$, $\sl \equiv 1$, ${\sigma_{n\ell}} \equiv 1$, and $f \equiv 1$. We observe that the domain configuration is such that $\Sigma$ intersects $\partial \Omega$ with a straight angle at the origin. Thus, according to Propositions \ref{['prop:local-regularity-mixed']} and \ref{['prop:regularity-small-s']}, in this setting, for $s < 1/2$, the regularity pickup in $\Omega_{1}$ is given by $\gamma(\Omega_{1}) = 3/2$. On the other hand, as $s \to 1$ one formally recovers a classical Poisson problem on the $L$-shaped $\Omega$ with homogeneous Dirichlet boundary conditions, that possesses a regularity pickup index $\gamma(\Omega) = 5/3$. We aim to illustrate the dependence of the solution regularity for this problem with respect to $s$. For this purpose, we used a mesh graded similarly to \ref{['eq:grading']}, although we considered the distance to the origin instead of the distance to the interface to drive the grading. The left panel in Figure \ref{['fig:slices_s01']} depicts the mesh in $\Omega$: we have a local mesh size near the origin $h \approx 1.5 \times 10^{-5}$. Singular behavior in $\Omega_{1}$ of solution for $s=0.1$. Left: top view; center: slice $x_1=0$; right: slice $x_2=0$. We performed least-squares fittings of our computed solution over the segments $\{0\} \times (0,1/200)$ and $(0,1/200) \times \{0\}$ with a model of the form $u(t) \simeq C t^\gamma$ using $11$ equally spaced sampling points, as shown in the center and right panels of Figure \ref{['fig:slices_s01']} for $s=0.1$; our findings are summarized in Table \ref{['tab:least-squares-graded']}. We observe that, for small $s$ and even up to $s \approx 0.7$, the asymptotic behavior of $u|_{\Omega_{1}}$ near the origin is consistent with the value $\gamma(\Omega_{1}) = 3/2$, there is a transition for larger values of $s$, and for $s\to 1$ one recovers an exponent consistent with the expected global regularity $\gamma(\Omega) = 5/3$. $s = 0.1$$s = 0.2$$s = 0.3$$s = 0.4$$s = 0.5$$s = 0.6$$x_1=0$$0.5046$$0.5041$$0.5037$$0.5032$$0.5028$$0.5028$$x_2=0$$0.5009$$0.5012$$0.5015$$0.5017$$0.5018$$0.5017$ $s = 0.7$$s = 0.8$$s = 0.9$$s=0.95$$s=0.99$$s=0.999$$x_1=0$$0.5046$$0.5142$$0.5544$$0.6007$$0.6527$$0.6655$$x_2=0$$0.5016$$0.5055$$0.5382$$0.5862$$0.6482$$0.6650$ Computed singular exponents via least-squares fitting on the segments $\{0\} \times (0,1/200)$ and $(0,1/200) \times \{0\}$ for different values of $s$ on a graded mesh. For the next experiment we add a negative, zero-order term to the energy $E_{II}$. More precisely, for $\omega > 0$ we consider $E_\omega (v) := E_{II}(v) - \omega^2 \| v \|_{L^2(\Omega)}^2$. Clearly, the energy minimization problem becomes ill-posed provided $\omega^2$ coincides with an eigenvalue of the associated differential operator over $\Omega$ with homogeneous boundary conditions. We test our code in the same setting as in §\ref{['sec:two-balls-experiment']}: since $\sl \equiv 1$ and ${\sigma_{n\ell}} \equiv 1$, formally, for $s \sim 1$ we expect to lose uniqueness of solutions as $\omega^2$ approaches any eigenvalue of the (classical) Laplacian with homogeneous Dirichlet boundary conditions. Over $\Omega = B(0,2) \subset {\mathbb R}^2$, such eigenvalues are given by $\lambda_{m,n} = \frac{j_{m,n}^2}{4},$ where $j_{m,n}$ is the $n$-th positive zero of the Bessel function of the first kind $J_m$. In increasing order, the first four numerical values are $\lambda_{0,1} \approx 1.45,\quad \lambda_{1,1} \approx 3.67, \quad \lambda_{2,1} \approx 6.59, \quad \lambda_{0,2} \approx 7.62.$ Figure \ref{['fig:conditioning-non-coercive']} reports the estimated $\ell_2$-condition numbers of the resulting system matrices on quasi-uniform meshes with $3442$ elements on $\Omega_{1}$ and 2249 on $\Omega_{2}$ --this corresponds to a mesh size $h \approx 5\times10^{-2}$-- as a function of $\omega^2$ and for different values of $s$ and $\omega^2$ spanning $(0,10)$. For $s=0.99$, we observe a very good agreement with the corresponding eigenvalues of the limit local problem while as $s$ decreases we observe that the spectrum tends to shift to the left. As $s \to 0$, the eigenvalues cluster around 1 from above, see the close-up for $\omega^2 \in (0,4)$ and $s=0.25$. In blue, computed condition numbers in the setting of §\ref{['sec:eigenvalues']} for $s=0.25$ (left), $s=0.75$ (center), and $s=0.95$ (right). In the right panel, the dashed line depicts the condition numbers for the corresponding discrete local problems. Our final experiment aims to compare the two energies \ref{['eq:energy-I']} and \ref{['eq:energy-III']}, and highlight the different fashion in which both recover the limiting local energy as $s \to 1$. On the square $\Omega = (-1/2,1/2)^2$, we consider $u_{loc}(x_1,x_2) = (1/4 - x_1^2) (1/4 - x_2^2)$, which solves $-\Delta u_{loc} = f$ with homogeneous Dirichlet boundary conditions and $f(x_1,x_2) = -2 (x_1^2 + x_2^2) + 1$. Considering $\Omega_{1} = \Omega \cap \{ x_1 < 0 \}$, Figure \ref{['fig:energy-comparison']} depicts the minimizers $u_I$ and $u_{II}$ we obtained resp. for the two energies $E_I$ and $E_{II}$ for $s = 0.999$, as well as the corresponding solution to the local problem. While visually both solutions to the coupled local/nonlocal models seem to approximate the desired solution, it seems the minimizer of $E_{II}$ is closer to it, while the coupling for the minimizer of \ref{['eq:energy-I']} is weaker. To make this statement more precise, Table \ref{['tab:energy-comparison']} compares the $L^2(\Omega)$-norms of such discrepancies for values of $s$ close to $1$. While for the energy $E_{II}$ the convergence seems to be linear with respect to $(1-s)$, the behavior for $E_I$ is less clear and the discrepancies are of a larger magnitude. Discrete energy minimizers for $s=0.999$ (left and center, corresponding to $E_I$ and $E_{II}$, respectively), and solution to the limiting local problem. We use the same scale for the three plots. $s = 0.9$$s = 0.99$$s = 0.999$$E_I$$1.6 \times 10^{-2}$$1.3 \times 10^{-2}$$3.2 \times 10^{-3}$$E_{II}$$7.1 \times 10^{-3}$$6.4 \times 10^{-4}$$6.4 \times 10^{-5}$ Discrepancy between solutions to \ref{['eq:Euler-Lagrange-E_I-discrete']} and \ref{['eq:weak-EIIv3-discrete']} for different values of $s$ and $u_{loc}$ in the $L^2(\Omega)$-norm. This appendix collects a few well-known results about Besov spaces that we use in the paper. While it is customary to introduce such spaces by real interpolation between integer-order Sobolev spaces, here we employ an equivalent definition through difference quotients. In particular, we restrict to spaces with differentiability order between 0 and 1. Given a function $v \colon {\mathbb R}^n \to {\mathbb R}$ and a vector $h \in {\mathbb R}^n \setminus \{0\}$, we denote by $\tau_hv(x) := v(x+h)$ the translation of $v$ with vector $h$. Given $\Lambda \subset {\mathbb R}^n$ and $\rho > 0$, we write $\Lambda_\rho := \{ x \in \Lambda \colon d(x, \partial\Lambda)<\rho \}.$ Let $D$ be a fixed ball centered at the origin, $\delta \in (0,1)$, $p,q \in [1,\infty]$. Then, $|v|_{B^\delta_{p,q}(\Lambda)} := (q\delta(1-\delta) \int_{D} \frac{\| \tau_h v - v \|_{L^p(\Lambda_{|h|})}^q}{|h|^{n+q\delta}}dh )^{1/q}$ for $p,q\in[1,\infty)$ while for $q=\infty$ we let $|v|_{B^\delta_{p,\infty}(\Lambda)} := \sup_{h \in D} \frac{\| \tau_h v - v \|_{L^p(\Lambda_{|h|})}}{|h|^{\delta}},$ Importantly, when $p = 2$ we have the algebraic and topological equivalence $H^\delta(\Lambda) = B^\delta_{2,2}(\Lambda)$ for all $\delta \in {\mathbb R}$, and actually the same equivalence between Besov and Sobolev spaces holds for an arbitrary $p \in [1,\infty]$ as long as the differentiability order $\delta$ is not an integer. We also make use of Besov spaces of functions that are supported on domains. For that purpose, we define $\dot{B}^\delta_{p,q}(\Lambda) := \{ v \in B^\delta_{p,q}({\mathbb R}^n) \colon \hbox{supp} v \subset \overline \Lambda \}.$ The Besov scale on bounded Lipschitz domains possesses some relevant orderings with respect to its defining parameters. While the inclusion $B^{\delta_1}_{p,q_1}(\Lambda) \subset B^{\delta_0}_{p,q_0}(\Lambda)\hbox{if} 0 < \delta_0 < \delta_1,1 \le p \le \infty ,1 \le q_0, q_1 \le \infty,$ is well-known (for example, Triebel10), we are interested in a quantitative estimate. The following result can be proven by repeating exactly the same steps as in BoLiNo23. Let $\Omega \subset {\mathbb R}^n$ be a bounded Lipschitz domain, $p,q \in [1,\infty)$, $\delta\in (0,1)$, and $\epsilon \in (0,1-\delta)$. Then, $B^{\delta+\epsilon}_{p,\infty}(\Omega) \subset B^{\delta}_{p,q} (\Omega)$ with $\|v\|_{ B^{\delta}_{p,q}(\Omega)} \lesssim \left( \frac{\delta}{\epsilon} \right)^{\frac{1}{q}} \|v\|_{B^{\delta+\epsilon}_{p,\infty}(\Omega)} \quad\forall \, v\in B^{\delta+\epsilon}_{p,\infty}(\Omega).$ Finally, we recall another well-known fact: that Besov spaces are increasing with respect to the parameter $q$. This follows from the characterization of their seminorms as discrete sums, $| v |_{B^\delta_{p,q}(\Lambda)} \simeq \left\lbrace \left( \sum_{k=0}^\infty | 2^{k\delta} \sup_{|h| \le 2^{-k}} \| \tau_h v - v \|_{L^p(\Lambda_{|h|})}|^q \right)^{1/q}\quad \hbox{if} q < \infty,\sup_{k \ge 0} \, 2^{k\delta} \sup_{|h| \le 2^{-k}} \| \tau_h v - v \|_{L^p(\Lambda_{|h|})}\quad \hbox{if} q = \infty , \right.$ and the fact that $\ell^q$ spaces are increasing with respect to $q$. We refer to DVLo for details. Let $\Omega \subset {\mathbb R}^n$ be a bounded Lipschitz domain, $p, \in [1,\infty)$, $1 \le q \le q' \le \infty$, and $\delta\in (0,1)$. Then, $B^{\delta}_{p,q}(\Omega) \subset B^{\delta}_{p,q'} (\Omega)$ with $\|v\|_{ B^{\delta}_{p,q'}(\Omega)} \lesssim \|v\|_{B^{\delta}_{p,q}(\Omega)} \quad\forall \, v\in B^{\delta}_{p,q}(\Omega).$ In this appendix, we prove Lemma \ref{['lemma:mapping-Laplacian']} about the mapping properties of the weighted fractional Laplacian \ref{['eq:def-wifl']}. We note that, because $0<\alpha < 1-2s$, according to Definition \ref{['eq:def-Besov']}, it suffices to show that $\| (\tau_h - I) (-\Delta_{\sigma_{n\ell}})^s u \|_{L^2({\mathbb R}^n)} \lesssim |h|^{\alpha} \| u \|_{H^{\alpha+2s}({\mathbb R}^n)}.$ For that purpose, given $x \in {\mathbb R}^n$ we split $(\tau_h - I) (-\Delta_{\sigma_{n\ell}})^s u (x) = (-\Delta_{\sigma_{n\ell}})^s u (x+h) - (-\Delta_{\sigma_{n\ell}})^s u (x)= \! \int_{{\mathbb R}^n} ( {\sigma_{n\ell}}(x+h, y+h) - {\sigma_{n\ell}}(x+h, x+h) ) \left( \frac{u(x+h) - u(y+h) - u(x) + u(y)}{|x-y|^{n+2s}} \right) dy+ \int_{{\mathbb R}^n} {\sigma_{n\ell}}(x+h, x+h) \left( \frac{u(x+h) - u(y+h) - u(x) + u(y)}{|x-y|^{n+2s}} \right) dy+ \int_{{\mathbb R}^n} \left( {\sigma_{n\ell}}(x+h, y+h) - {\sigma_{n\ell}}(x, y) \right) \left( \frac{u(x) - u(y)}{|x-y|^{n+2s}} \right) dy=: I(x) + II(x) + III(x) .$ We begin bounding the term $I$. We exploit the boundedness and Hölder regularity of ${\sigma_{n\ell}}$, cf. \ref{['eq:Holder-sigma']}, and fix some $\theta < \min\{2s, \beta\}$, to deduce $|{\sigma_{n\ell}}(x+h, y+h) - {\sigma_{n\ell}}(x+h, x+h) | \lesssim |x-y|^\theta.$ With that, making the change of variables $z = y - x$ and using Minkowski's inequality, $\| I \|_{L^2({\mathbb R}^n)}\lesssim \left( \int_{{\mathbb R}^n} \left(\int_{{\mathbb R}^n} \frac{|u(x+h)-u(y+h) -u(x) + u(y)|}{|x-y|^{n+2s-\theta}} \, dy \right)^2 dx \right)^{1/2}= \left( \int_{{\mathbb R}^n} \left(\int_{{\mathbb R}^n} \frac{|u(x+h)-u(x+h+z) -u(x) + u(x+z)|}{|z|^{n+2s-\theta}}\, dz \right)^2 dx \right)^{1/2}\le \int_{{\mathbb R}^n} \frac{\| \tau_z( u - \tau_h u) - ( u - \tau_h u)\|_{L^2({\mathbb R}^n)}}{|z|^{n+2s-\theta}} \, dz\lesssim \| u - \tau_h u \|_{B^{2s-\theta}_{2,1}({\mathbb R}^n)}.$ Next, we observe that, by interpolation of the inequalities $\| u - \tau_h u \|_{L^2({\mathbb R}^n)} \lesssim |h|^{\alpha+2s} \| u \|_{H^{\alpha+2s}({\mathbb R}^n)} , \qquad \| u - \tau_h u \|_{H^{\alpha+2s}({\mathbb R}^n)} \lesssim \| u \|_{H^{\alpha+2s}({\mathbb R}^n)} ,$ we have $\| u - \tau_h u \|_{B^{2s-\theta}_{2,1}({\mathbb R}^n)} \lesssim |h|^{\alpha + \theta} \| u \|_{H^{\alpha+2s}({\mathbb R}^n)}.$ Therefore, since $|h| \le 1$, we obtain $\| I \|_{L^2({\mathbb R}^n)} \lesssim |h|^{\alpha + \theta} \| u \|_{H^{\alpha+2s}({\mathbb R}^n)} \le |h|^{\alpha} \| u \|_{H^{\alpha+2s}({\mathbb R}^n)} .$ We now deal with the term $II$ in \ref{['eq:split-translation-Laplacian']}. For this, we note that it essentially coincides with a first difference for a non-weighted fractional Laplacian, namely, $II (x) = {\sigma_{n\ell}}(x+h, x+h) \left( \tau_h - I \right) (-\Delta)^s u (x),$ and thus $\| II \|_{L^2({\mathbb R}^n)} \lesssim \| \left( \tau_h - I \right) (-\Delta)^s u \|_{L^2({\mathbb R}^n)} \lesssim |h|^{\alpha} \| (-\Delta)^s u \|_{B^{\alpha}_{2,\infty}({\mathbb R}^n)}.$ We use the continuity of the embedding $H^{\alpha}({\mathbb R}^n) \subset B^{\alpha}_{2,\infty}({\mathbb R}^n)$ (see \ref{['eq:Besov-Besov-emb']}) and the well-known mapping property $(-\Delta)^s \colon H^{\alpha+2s}({\mathbb R}^n) \to H^{\alpha}({\mathbb R}^n)$ to obtain $\| II \|_{L^2({\mathbb R}^n)} \lesssim |h|^{\alpha} \| u \|_{H^{\alpha+2s}({\mathbb R}^n)}.$ Finally, we bound the integral $III$ in \ref{['eq:split-translation-Laplacian']}. We again exploit the $\beta$-Hölder regularity of ${\sigma_{n\ell}}$ \ref{['eq:Holder-sigma']}, $\| III \|_{L^2({\mathbb R}^n)} \lesssim |h|^\beta \left( \int_{{\mathbb R}^n} \left( \int_{{\mathbb R}^n} \frac{|u(x) - u(y)|}{|x-y|^{n+2s}} \, dy \right)^2 dx \right)^{1/2} .$ Upon making the change of variables $z = x + y$ and using Minkowski's inequality, we obtain $\left( \int_{{\mathbb R}^n} \left( \int_{{\mathbb R}^n} \frac{|u(x) - u(y)|}{|x-y|^{n+2s}} \, dy \right)^2 dx \right)^{1/2} = \int_{{\mathbb R}^n} \frac{\| \tau_z u - u \|_{L^2({\mathbb R}^n)}}{|z|^{n+2s}} \, dz$ and therefore, by definition \ref{['eq:def-Besov']}, $\| III \|_{L^2({\mathbb R}^n)} \lesssim |h|^{\beta} \| u \|_{B^{2s}_{2,1}({\mathbb R}^n)}.$ The bound \ref{['eq:mapping-Laplacian']} follows by putting together \ref{['eq:bound-I']}, \ref{['eq:bound-II']}, and \ref{['eq:bound-III']} and noticing that $\| u \|_{B^{2s}_{2,1}({\mathbb R}^n)} \lesssim \| u \|_{H^{\alpha+2s}({\mathbb R}^n)}$ by \ref{['eq:Besov-Sobolev-emb']}-\ref{['eq:Besov-Besov-emb']}. G. Acosta, F. Bersetche, and J. Borthagaray. A short FE implementation for a 2d homogeneous Dirichlet problem of a fractional Laplacian. Comput. Math. Appl., 74(4):784--816, 2017.G. Acosta, F. Bersetche, and J. D. Rossi. Local and nonlocal energy-based coupling models. SIAM Journal on Mathematical Analysis, 54(6):6288--6322, 2022.G. Acosta, F. Bersetche, and J. D. Rossi. Coupling local and nonlocal equations with Neumann boundary conditions. Revista de la Unión Matemática Argentina, 65(2):533--565, 2023.G. Acosta and J. P. Borthagaray. A fractional Laplace equation: Regularity of solutions and finite element approximations. SIAM J. Numer. Anal., 55(2):472--495, 2017.D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert. Application of a fractional advection-dispersion equation. Water Resources Research, 36(6):1403--1412, 2000.K. Bogdan, T. Grzywny, K. Pietruska-Pałuba, and A. Rutkowski. Extension and trace for nonlocal operators. J. Math. Pures Appl. (9), 137:33--69, 2020.J. P. Borthagaray and P. Ciarlet Jr. Nonlocal models for interface problems between dielectrics and metamaterials. In International Congress on Engineered Material Platforms for Novel Wave Phenomena, 2017.J. P. Borthagaray, W. Li, and R. H. Nochetto. Quasi-linear fractional-order operators in Lipschitz domains. SIAM Journal on Mathematical Analysis, 56(3):4006--4039, 2024.J. P. Borthagaray and R. H. Nochetto. Besov regularity for the Dirichlet integral fractional Laplacian in Lipschitz domains. Journal of Functional Analysis, 284(6):109829, 2023.J. P. Borthagaray and R. H. Nochetto. Constructive approximation on graded meshes for the integral fractional Laplacian. Constructive Approximation, 57(2):463--487, 2023.G. Capodaglio, M. D'Elia, P. Bochev, and M. Gunzburger. An energy-based coupling approach to nonlocal interface problems. Computers & Fluids, 207:104593, 2020.Z. Chen and R. H. Nochetto. Residual type a posteriori error estimates for elliptic obstacle problems. Numerische Mathematik, 84:527--548, 2000.P. Ciarlet, Jr. Analysis of the Scott-Zhang interpolation in the fractional order Sobolev spaces. J. Numer. Math., 21(3):173--180, 2013.P. Clément. Approximation by finite element functions using local regularization. Rev. Française Automat. Informat. Recherche Opérationnelle Sér. RAIRO Analyse Numérique, 9(R-2):77--84, 1975.R. Cont and P. Tankov. Financial modelling with jump processes. Chapman & Hall/CRC Financial Mathematics Series. Chapman & Hall/CRC, Boca Raton, FL, 2004.M. Costabel, M. Dauge, and S. Nicaise. Corner singularities and analytic regularity for linear elliptic systems. part I: Smooth domains. 2010.M. Dauge. Neumann and mixed problems on curvilinear polyhedra. Integral equations and operator theory, 15:227--261, 1992.R. A. DeVore and G. G. Lorentz. Constructive approximation, volume 303. Springer Science & Business Media, 1993.P. Diehl and S. Prudhomme. Coupling approaches for classical linear elasticity and bond-based peridynamic models. Journal of Peridynamics and Nonlocal Modeling, 4(3):336--366, 2022.S. Dipierro, X. Ros-Oton, and E. Valdinoci. Nonlocal problems with Neumann boundary conditions. Revista Matemática Iberoamericana, 33(2):377--416, 2017.K. Disser, H.-C. Kaiser, and J. Rehberg. Optimal sobolev regularity for linear second-order divergence elliptic operators occurring in real-world problems. SIAM Journal on Mathematical Analysis, 47(3):1719--1746, 2015.Q. Du, X. H. Li, J. Lu, and X. Tian. A quasi-nonlocal coupling method for nonlocal and local diffusion models. SIAM Journal on Numerical Analysis, 56(3):1386--1404, 2018.M. D’Elia, X. Li, P. Seleson, X. Tian, and Y. Yu. A review of local-to-nonlocal coupling methods in nonlocal diffusion and nonlocal mechanics. Journal of Peridynamics and Nonlocal Modeling, pages 1--50, 2022.C. Ebmeyer. Mixed boundary value problems for nonlinear elliptic systems in n-dimensional Lipschitzian domains. Zeitschrift für Analysis und ihre Anwendungen, 18(3):539--555, 1999.C. Ebmeyer and J. Frehse. Mixed boundary value problems for nonlinear elliptic equations in multidimensional non-smooth domains. Mathematische Nachrichten, 203(1):47--74, 1999.G. Estrada-Rodriguez, H. Gimperlein, and J. Stocek. Nonlocal interface problems: Modeling, regularity, finite element approximation. Preprint, http://www.macs.hw.ac.uk/hg94/interface.pdf, 2020.B. Faermann. Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary element methods. I. The two-dimensional case. IMA J. Numer. Anal., 20(2):203--234, 2000.B. Faermann. Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary element methods. II. The three-dimensional case. Numer. Math., 92(3):467--499, 2002.C. G. Gal and M. Warma. Nonlocal transmission problems with fractional diffusion and boundary conditions on non-smooth interfaces. Communications in Partial Differential Equations, 42(4):579--625, 2017.H. Gimperlein, E. Stefan, and J. Stocek. Corner singularities for the fractional laplacian and finite element approximation. Preprint, http://www.macs.hw.ac.uk/hg94/corners.pdf, 2019.P. Grisvard. Behavior of the solutions of an elliptic boundary value problem in a polygonal or polyhedral domain. In Numerical Solution of Partial Differential Equations--III, pages 207--274. Elsevier, 1976.P. Grisvard. Elliptic problems in nonsmooth domains, volume 69 of Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, 2011.F. Grube and T. Hensiek. Robust nonlocal trace spaces and Neumann problems. Nonlinear Analysis, 241:113481, 2024.F. Grube and M. Kassmann. Robust nonlocal trace and extension theorems. arXiv preprint arXiv:2305.05735, 2023.Q.-Y. Guan. Integration by parts formula for regional fractional Laplacian. Communications in Mathematical Physics, 266(2):289--329, 2006.D. Jerison and C. E. Kenig. The inhomogeneous Dirichlet problem in Lipschitz domains. Journal of functional analysis, 130(1):161--219, 1995.F. Jochmann. An $H^s$-regularity result for the gradient of solutions to elliptic equations with mixed boundary conditions. Journal of mathematical analysis and applications, 238(2):429--450, 1999.M. Kassmann and W. Madych. Difference quotients and elliptic mixed boundary value problems of second order. Indiana University mathematics journal, pages 1047--1082, 2007.V. A. Kondrat'ev. Boundary value problems for elliptic equations in domains with conical or angular points. Trudy Moskovskogo Matematicheskogo Obshchestva, 16:209--292, 1967.D. Kriventsov. A local-nonlocal transmission problem. PhD thesis, University of Texas at Austin, 2015.D. Kriventsov. Regularity for a local--nonlocal transmission problem. Archive for Rational Mechanics and Analysis, 217(3):1103--1195, 2015.M. Kwaśnicki. Ten equivalent definitions of the fractional Laplace operator. Fract. Calc. Appl. Anal., 20(1):7--51, 2017.P.-A. Raviart and J. M. Thomas. Introduction à l'analyse numérique des équations aux dérivées partielles. Masson, 1983.X. Ros-Oton and J. Serra. The Dirichlet problem for the fractional Laplacian: regularity up to the boundary. J. Math. Pures Appl. (9), 101(3):275--302, 2014.G. Savaré. Regularity results for elliptic equations in Lipschitz domains. J. Funct. Anal., 152(1):176--201, 1998.F. J. Sayas, T. S. Brown, and M. E. Hassell. Variational techniques for elliptic partial differential equations: Theoretical tools and advanced applications. CRC Press, 2019.C. Schneider and F. O. Szemenyei. Sobolev meets Besov: Regularity for the Poisson equation with Dirichlet, Neumann and mixed boundary values. Analysis and Applications, 20(05):989--1023, 2022.L. R. Scott and S. Zhang. Finite element interpolation of nonsmooth functions satisfying boundary conditions. Math. Comp., 54(190):483--493, 1990.S. A. Silling. Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys. Solids, 48(1):175--209, 2000.H. Triebel. Theory of Function Spaces. Modern Birkhäuser Classics. Springer Basel, 2010.M. Warma. The fractional relative capacity and the fractional Laplacian with Neumann and Robin boundary conditions on open sets. Potential Anal., 42(2):499--547, 2015.