Singular boundary behaviour and large solutions for fractional elliptic equations

Abstract We perform a unified analysis for the boundary behaviour of solutions to nonlocal fractional equations posed in bounded domains. Based on previous findings for some models of the fractional Laplacian operator, we show how it strongly differs from the boundary behaviour of solutions to elliptic problems modelled upon the Laplace–Poisson equation with zero boundary data. In the classical case it is known that, at least in a suitable weak sense, solutions of the homogeneous Dirichlet problem with a forcing term tend to zero at the boundary. Limits of these solutions then produce solutions of some non‐homogeneous Dirichlet problem as the interior data concentrate suitably to the boundary. Here, we show that, for equations driven by a wide class of nonlocal fractional operators, different blow‐up phenomena may occur at the boundary of the domain. We describe such explosive behaviours and obtain precise quantitative estimates depending on simple parameters of the nonlocal operators. Our unifying technique is based on a careful study of the inverse operator in terms of the corresponding Green function.


INTRODUCTION
In recent years, there have been many studies on boundary value problems driven by nonlocal operators L obtained as fractional powers of uniformly elliptic operators, such as the Laplacian. In this context, according to the 'degree of nonlocality' of the leading operator in the differential equation, additional values need to be prescribed either on the boundary of the underlying domain or on its whole complement. So, given a regular bounded domain Ω ⊆ ℝ , the simplest Dirichlet problems take the form of an equation the last choice depending on the nonlocal operator L. Sometimes, (1.1) is written in some weak form which also encodes (1.2). In the standard elliptic theory, (1.2) can be replaced by = g on Ω, for g an function and therefore a.e. finite on Ω. In this paper, we study solutions to equations of the form (1.1) that develop an explosive behaviour at the boundary, that is, solutions satisfying ( ) → +∞, as → 0 , for almost all 0 ∈ Ω.
They are usually called large solutions and they account for a new phenomenon, not appearing in the classical elliptic theory. We will show that large solutions are tightly connected to the solutions of the homogeneous problem via a natural limiting process. Finally, they exhibit quite peculiar divergence rates that we will derive. All of this will be done for a specific class of nonlocal operators L that includes the usual examples and more. The present research is motivated by two striking results involving singular behaviour near the boundary for the solutions of (1.1) in the case where L is the so-called restricted fractional Laplacian (RFL; for which one has to prescribe data on ℝ ⧵ Ω).
One of these striking results is the existence of nontrivial solutions of (1.1) such that = 0 in Ω which, moreover, are positive everywhere and blow up on the boundary. Explicit examples on the ball were constructed in [33] (see also [6,8]). The existence of this kind of solutions was systematised independently in [30] (which also contains a thorough regularity theory, see also [29] for related results and [31] for a review) and in [1] for the RFL, and extended in [2] for the spectral fractional Laplacian (SFL; which requires prescribed data at the boundary). These correspond to positive harmonic functions in the theory of the standard Laplacian, although the classical theory does not admit any harmonic function with uniform blow-up at the boundary.
The second striking result, described in [1] when L is the RFL and in [2] when L is the SFL, is that some admissible functions produce solutions blowing-up at the boundary, although they are limits of solutions with 'nice' and zero boundary data. This is as well a new behaviour of the nonlocal problem, not present for the usual Laplacian.
For the case of the usual Laplacian, it is known that the wider classes of weak or very weak solutions obtained as limits of the variational solutions satisfy the boundary condition either in the sense of traces or in a more generalised sense, described in [36] as the average condition The aim of the present work is to show that these two blow-up phenomena occur for a large class of nonlocal operators of elliptic type. We treat in a unified way the typical nonlocal elliptic equations, in particular the different fractional Laplacians on bounded domains. Our distinctive technique is based on the use of the Green kernel which gives a common roof to the several different cases. This approach extends previous work in [10,28].
We consider a general family of operators indexed on two parameters: one describing the interior point singularity of the Green kernel, the other one the kernel's boundary behaviour. This requires serious technical work, that justifies the extension of the paper.
First, we want to study and classify the explosive (or large) solutions whose singularity is, in some sense, generated by the right-hand side . In particular, we compute explicitly the asymptotic boundary behaviour of for the family of power-like data ≍ near the boundary, where ( ) ∶= dist( , Ω). Here, we say that ≍ g on a set if there exists > 0 such that −1 g ⩽ ⩽ g on that set. Our main formula (1.5) gives the behaviour of in terms of and the kernel of L in simple algebraic terms. The formula covers the whole range of behaviours, explosive or not. We also translate estimate (1.3) to our context by introducing a suitable weight, taking care of the singular profiles (see Lemma 3.8). We provide some careful numerical computations, to show the formation of the boundary singularity due to the right-hand side (see Figures 3 and 4).
Even if the solution operator for the Dirichlet problem acting on a class of good functions produces solutions with Dirichlet boundary data, we show that the natural closure of that solution operator to its maximal domain of definition produces solutions which no longer satisfy the Dirichlet condition and could reach a range of boundary blow-up that we describe. In the case of problems in which the boundary condition is set (for example, the SFL), this is counterintuitive. The occurrence of boundary blow-up is a very important fact, that does not happen for the usual Laplacian.
Second, we remark how there is a different class of explosive solutions whose singularity is not generated by any right-hand side. In fact, they can be chosen as 'L-harmonic in Ω' in the sense L = 0. This class relies on some hidden information in the form of singular behaviour that can be prescribed on the boundary. Moreover, this second class can be obtained as a limit of singular solutions of the previous class as the support of concentrates at the boundary in a convenient way. This means they cannot be disregarded in any complete theory of the problem. See the detailed results in Section 4.
We conclude this introduction with an important remark. If a definition of solution of (1.1) is 'too weak', then the combination of the two classes seems to pose a problem to uniqueness, as it happens in the classical elliptic theory. This highlights the importance of a suitable definition of weak solution of (1.1) preserving uniqueness and including the classical solutions. We provide this definition in Section 2 under the name of weak-dual solution, and show that the problem is then well-posed. We also detect the optimal class of admissible data . To take care of the second class, we construct a 'singular boundary data' problem. We give a well-posed notion of solution for this second problem: uniqueness is the easy part.

Main topics and results
Existence and uniqueness results for (1.1) To begin with, we need to produce a general existence theory for data in good classes, that is, compactly supported and bounded. We want to treat a general class of operators such that the unique solution of (1.1)-(1.2) is given by the formula with kernels ∶ Ω × Ω → ℝ such that, for any , ∈ Ω, ( , ) = ( , ), Ω ∈ 1,1 , and (K1) .
The two exponents and take values Relation of parameters and such that ( ) ≍ Their relative values will play an important role in the results. Note that hypothesis (K1) means that  is self-adjoint. This allows to cover a large class of self-adjoint operators L, as we will describe in Subsection 1.2. Throughout this note, we use the notation ∧ = min{ , }, ∨ = max{ , }.
The main structural assumption on is (K2). This general assumption was introduced in [11] to cover some notable examples of this general class of operators given by the three most known fractional Laplacian operators: (i) The RFL: in this case = ∈ (0, 1).
These examples will be presented in some more detail in Subsection 1.2, so that we can adapt to them the general results. As we will explain, the estimate (K2) in each of these examples is recovered by ad hoc techniques in different papers. In Section 5, we will make a brief discussion on the general setting.
In Section 2, we prove existence, uniqueness, a priori estimates, and some regularity for problem (1.1). In Section 3, we prove that the optimal class of data such that (K0) is well-defined (meaning (| |) ≢ +∞) is

Boundary behaviour
As we mentioned, for the standard Laplacian −Δ, the zero boundary data are taken in some sense even when is taken in the optimal class of data. The sense depends on how good is , see the general results in [36]. A quite novel property of the RFL on bounded domains shows that this is not true for admissible even if they are not so badly behaved. This is explained in [1] and we want to extend the analysis to our general class of operators and show the detailed relation between the operators, the boundary behaviour of , and the singular boundary behaviour of the solution. The main information about the operators will be the values of and . In Theorem 3.4, we establish the explicit estimate ( ) ≍ ∧( +2 ) whenever + > −1 and ≠ − 2 (1.5) that needs a delicate computation using the properties of the kernel. This is depicted in Figure 1. Note that + > −1 is the condition, so that = belongs to the admissible class given by (1.4).

F I G U R E 2 Different relations between and
In many cases, the existence of eigenfunctions is known, and their boundary behaviour is wellunderstood. Under (K0), (K2), and some extra assumptions on the operator L, the authors in [10] proved that the operator  admits an eigendecomposition and its first eigenfunction Φ 1 satisfies The boundary behaviour is clear from the algebraic point of view, since is the only exponent fixed by .

Solutions with singular behaviour
We observe that, according to formula (1.5), there are values of for which the solution associated to datum is singular at the boundary: this happens whenever ∈ (−1 − , −2 ) is allowed, and therefore when > 2 − 1. In particular, it comes out that if > 2 − 1 , then there exist solutions of the Dirichlet problem not complying with the condition = 0 on the boundary. This was known for the RFL [1, Proposition 3] and the SFL [23,Proposition 7], whereas this phenomenon does not take place for the CFL (as we show below). The behaviour ∧(2 − −1) , corresponding to the limit case = −1 − , serves somehow as an upper bound for solutions. In Lemma 3.8, we will prove that (a) if > − 1∕2, then for any ∈ 1 (Ω, ) (b) if < − 1∕2, then for any nonnegative ∈ 1 (Ω, ) and > 0 We also prove that, in the case = − 1∕2, there is a logarithmic correction. For the usual Laplacian, when = = 1, we have 0 = 2 − − 1 < : this reproduces (1.3). This same fact holds for the CFL, because = 2 − 1. If 2 − − 1 > 0, then all solutions tend to 0 upon approaching the boundary.
The two conditions > 2 − 1 and > − 1 2 allow us to split the parameter , in three regions as in Figure 2.

Normal derivatives
A sharper study of the boundary behaviour of solution with data ∈ ∞ (Ω) consists of the analysis of the limit We will call this limit -normal derivative. We devote Subsection 3.3 to the study of these normal derivatives (see Theorem 3.15).

Large solutions
In [8], the authors introduce a surprising singular solution of the homogeneous problem = 0 that shows very precise asymptotics at the boundary. It is the type known as large solution in other situations for nonlinear equations. For example, the function is known to satisfy (−Δ) RFL ( ) = 0 for | | < 1, see [6, Example 1] and [33]. In [1], there is a complete description of the singular boundary value problem for the RFL, while in [2] there is the analogue for the SFL. Note that, in the limit → 1, the in (1.6) becomes the characteristic function of the ball. We prove that this theory may be obtained as a limit of interior problems. We construct one such particular large solution ⋆ which is L-harmonic on the interior (L ⋆ = 0 in Ω). In Section 4, we show that there exists a sequence of admissible functions ( ) ∈ℕ (with dist(supp , Ω) < 2∕ ) such that This limit function has the boundary behaviour except in the case = − 1∕2 when a logarithmic correction is in order. Note that the exponent is the upper bound of the range in (1.5). We will prove that the problems † for any ∈ ∞ (Ω).
have a unique solution , which is comparable to ⋆ at the boundary. Going back to the representation in L, this would mean that = [ ] and then formally Note that this is the formulation of the Dirichlet boundary value problem when L = −Δ so = = 1. When > − 1 2 , this is even the unique (see Theorems 4.6 and 4.13) weak-dual solution of (1.8)

Comments
Our presentation unifies in a single theory previous results for the RFL ( ⋆ ≍ −1 , see [1]), the SFL ( ⋆ ≍ 2( −1) , see [2]), the CFL ( ⋆ ≍ 1, see [14]), and even the classical Laplacian ( ⋆ ≍ 1). The case < − 1 2 , which does not include any of the main known examples, is somewhat particular. In this case, due to (1.7), 0 < ⋆ ( ) → 0 as → Ω and it is a non-trivial solution of (1.1) with data = 0. This yields some doubt about the uniqueness of solutions to (1.1)-(1.2). Furthermore, if ⩽ − 1 2 , then ⋆ ≍ , which in turn means that the critical solutions have the same boundary behaviour as the solutions for regular data . This does not seem to be consistent with elliptic problems like (1.1).

Some examples
Large classes of operators L have Green operators  given by (K0)-(K3): here are some notorious examples that are reviewed, for instance, in [3,10,40].

The restricted fractional Laplacian
The RFL is defined as the singular integral operator up to a multiplicative constant only depending on and , and corresponds to the power of the Laplacian operator defined in ℝ (which can be equivalently defined via the Fourier transform). The natural boundary conditions are given in ℝ ⧵ Ω and also have a stochastic interpretation corresponding to killing a Lévy flight upon leaving Ω.
Here we can consider all ∈ (0, 1) and we have the precise value = .
Details can be consulted in many references, see, for instance, [11,37].
Perturbations of the RFL Using the above one, it is possible to build other examples. Here are a couple of interesting operators which are included in our analysis and the corresponding references: • (−Δ) RFL + ⋅ ∇ for ∈ (1∕2, 1) and ∈ ∞ (Ω): in this case (see [9]); = .

The spectral fractional Laplacian
A different way of considering the power of the Laplacian consists in taking the power of the Dirichlet Laplacian, that is, the Laplacian coupled with homogeneous boundary conditions. This approach typically makes use of an eigenbasis expansion. Let ( ) ∈ℕ be the eigenfunctions of the Laplacian linked to the nondecreasing sequence of eigenvalues 0 < 1 < 2 ⩽ … (repeated according to their multiplicity) Let ∈ 2 ∩ 1 0 (Ω). Lettingˆ= ∫ Ω , we have the representation The SFL is the operator with eigenvalues corresponding to eigenfunctions . Hence, we define Since this is an operator-wise definition we provide the boundary conditions given from the classical operator, and hence the problem is We underline how this is not the only possible representation and it is also possible to write it as the singular integral operator and Ω the Dirichlet heat kernel on Ω. It is possible to prove that, when Ω ∈ 1,1 , Stochastically speaking, this operator generates a subordinate killed Brownian motion, which is a Brownian motion killed upon hitting Ω and which is then evaluated at random times distributed as an increasing -stable process in (0, ∞), see [39]. The killing of the Brownian motion as it touches the boundary is encoded in the homogeneous boundary conditions.
Here again ∈ (0, 1) and in this case, Details can be consulted in many references, see, for instance, [11,13].
An interpolation of the RFL and the SFL A family of 'intermediate' operators between the RFL and the SFL has been built in [35]. For 1 , 2 ∈ (0, 1], one can consider the spectral 2 power of a RFL of exponent 1 It is formally clear that The Green function associated to this operator satisfies (K2) with

The censored fractional Laplacian
This operator is defined as (recall that the RFL is evaluated only on functions satisfying = 0 in ℝ ⧵ Ω). This operator generates a censored stable process, introduced in [7], a stable process which is confined in Ω and finally killed upon hitting Ω. For this reason, a suitable boundary condition is Here ∈ (1∕2, 1) and see [7,15]. A class of operators which generalises and includes the CFL is covered by the analysis in [18].
Proof. We are going to extensively use assumptions (K0)-(K3) without further notice.
To prove (2.4), we note that, for any ⋐ Ω and ∈ 1 (Ω, ), we have and this proves the result. □ Remark 2.2. In Subsection 3.3, we will give a sharper characterisation of the image of map  in terms of weighted 1 spaces.
3. Formally, one could take ∈ (Ω) and estimate where we have used (2.1) on ( Ω ) or take ∈ (Ω, ) and estimate, for any ⋐ Ω where we have used (2.2) on ( ). This computation is justified in the typical examples where is continuous. However, since we have made no continuity assumptions for , it is possible that integration against a measure is not defined. We will give more details on this case in Subsection 2.5.

Weak-dual formulation
If L is self-adjoint, equations of type (1.1) are typically written in very weak form as for all test functions in some adequate space given by the operator and the boundary conditions. Since we want to tackle multiple types of operators and boundary conditions, we focus instead on the weak-dual formulation (see, for example, [10]). This is formulated instead in terms of the inverse operator , which is taken as an a priori. This allows to avoid giving a meaning to L .
Note that this weak-dual formulation is equivalent to take test functions ∈ ( ∞ (Ω)) in (2.8). Also, we underline how the integral in the right-hand side of (2.9) is finite in view of (2.2).
This function is precisely = ( ) and it satisfies for any ⋐ Ω.
Proof. Let us first note that = ( ) ∈ 1 loc (Ω) in view of (2.4). It formally satisfies (2.10) as a consequence of (K1) by the Fubini's theorem. This formal bounds are indeed rigourous for ∈ ∞ (Ω). Furthermore, due to the bounds provided by Theorem 2.1 one can pass to the limit in approximations.
We now focus on uniqueness. Let 1 , 2 be two solutions to (2.10). Then Let ⋐ Ω and = sign( 1 − 2 ) ∈ ∞ (Ω). Using this as a test function, we deduce Since this holds for every ⋐ Ω, we have that 1 = 2 a.e. in Ω. Also, we have that which is a nontrivial inequality thanks to (2.2). □
Thanks to Theorem 2.6 and Remark 2.8, we have shown that 1 (Ω, ) is the optimal class of data.

Uniform integrability over compacts
Let us show that  maps 1 -bounded sequences into 1 -weakly pre-compact sequences.
In particular, for any ⋐ Ω,  maps bounded sequences in 1 (Ω, ) into uniformly integrable sequences in .
Proof. We have that We take 1 < < −2 . Due to the Hölder's inequality, We estimate this last integral to recover . For any ∈ Ω such that dist( , ) < dist( , Ω)∕2 we have that where depends only on and Ω. One the other hand, if is such that dist( , ) ⩾ dist( , Ω)∕2, for ∈ we have | − | > dist( , ) and so we compute where depends only on and Ω. This completes the proof. □

Measure data and continuous solutions
Under mild assumptions on the Green kernel , it is possible to improve (2.1) and (2.2) to higher regularity of solutions. By duality, this allows more general data in (1.1) and we are particularly interested in measure data. For this reason, let us assume that Proof. In view of (2.1) and (2.2), we just need to justify the continuity claim. Let us consider ∈ ∞ (Ω). To prove continuity we select an ∈ Ω, and ( ) ∈ℕ ⊂ Ω such that → as ↑ ∞. By assumption (K4) we know that ( , ⋅) → ( , ⋅) a.e. in Ω. Moreover, let us note and Ω is bounded. Therefore, so is ( ( , ⋅) ) ∈ℕ . Due to the weak compactness in reflexive spaces it is convergent. Applying (K4) we can compute the pointwise limit This proves that ( ) ∈ (Ω). Let ∈ ∞ (Ω). Since we have already proven that ( ) ∈ (Ω), we only need to prove that − ( ) is continuous on some small neighbourhood of Ω. Consider > 0 small enough Select now an ∈ and let Ω ∈ → as ↑ +∞. Since is open, then ∈ for large enough. Again, by weak compactness, This completes the proof. □ With this new machinery, we can justify the intuition given by Remark 2.3.
Let ⋐ Ω and = sign( 1 − 2 ) ∈ ∞ (Ω). Using this as a test function, we deduce Since this holds for every ⋐ Ω, we have that 1 = 2 a.e. in Ω. Also, we have that which is a nontrivial inequality thanks to (2.2). □

BREAKDOWN OF THE BOUNDARY CONDITION IN THE INTERIOR PROBLEM
We address now the main question of this paper, which is the violation of the boundary data in the optimal theory for the interior problem. We give precise answers of the anomalous boundary behaviour in terms of the behaviour of the forcing data.

Range of exponents
Before stating and proving the main result of this paragraph, we need to state a couple of technical estimates on which the result is based. Since the proofs of these estimates is rather long and technical, we defer them to Appendix B. The first one gives some interior estimates; the second one is describing the sharp behaviour of solutions at the boundary.
In what follows, we use without further note > 0 to denote the fixed width on which the tubular neighbourhood theorem can be rightfully applied, that is, the map defines a diffeomorphism to its image. Here, represents the interior unit normal. This is well-known for smooth manifolds (see [21]), and holds also for 1,1 open sets of ℝ . The notation might seem like an abuse of notation, but it will lead to no confusion since, in this setting, dist(Φ( , ), Ω) = for sufficiently small.  For the RFL we take the weights for the finite difference discretisation of the fractional Laplacian in ℝ in [22] (see also [19]). The discretisation for smooth functions is rigorously shown to be (ℎ 2 ). A previous approach by Finite Differences is given in [34]. Experimental results in [34] suggest that the restriction for RFL is of order (ℎ ). For the SFL we use as a discretisation the fractional power of the finite differences matrix of the usual Laplacian (−Δ): it is known that the eigenvalues of this matrix converge to those of the usual Laplacian, and hence its fractional power produces a convergent scheme for the SFL. A different scheme can be found in [20]. For the CFL we have used a novel approach, which we will describe in an upcoming paper. (3.6) We are now ready to prove the following estimate.
Equation (3.7) can be interpreted by means of formula (1.5) and Figure 1.
Proof. Let us first note that conditions ⩽ − 2 and ⩽ − 1 2 are not compatible: indeed, if they both held, then it would be + ⩽ 2 − 2 ⩽ −1 contradicting our standing assumption on .
We pick some fixed < and we write Note that For the other term, we exploit Lemma 3.3 and (3.1) to say where Θ is defined as in Lemma 3.3. Now, the asymptotic behaviour is driven by the least exponent on , yielding the situation depicted in (3.7). □ Remark 3.5. Let us look at the ranges for and as in Theorem 3.4, disregarding the logarithmic cases, to better understand the possible boundary behaviours of solutions to (1.1). When > − 1 2 the admissible range for is (−1 − , +∞); in this case runs in (2 − − 1, ]: note that 2 − − 1 might be negative, meaning that also is allowed to be negative in some cases. This translates in particular into a rebuttal of ( ) = 0 on Ω, despite the fact that this would be the solution to a homogeneous boundary (or exterior problem) value problem. For exterior problem, this shows solutions are discontinuous on the boundary for some singular data (possibly outside 1 (Ω)). For boundary value problems, this is a breakdown of the boundary condition. However, this behaviour intrinsic to the problem, since we are only constructing the closure of the solution operator , to its maximal domain of definition. If instead < − 1 2 , then again ranges in (− − 1, +∞), but this time is bound to be equal to , meaning that there is no range for . Example 3.6. Let us exemplify the statement of Theorem 3.4. If we consider = 0, we deduce in Ω.

Below this value, if is of the form
Remark 3.7. Note that, if ∈ (−1∕2, −2 ) we have that ∈ 2 (Ω) and ( ) ∉ ∞ (Ω). This is possible if ∈ (0, 1∕4). Hence, this breakdown of the boundary conditions happens inside the variational (energy) theory. This should not be surprising since, for < 1∕2, 0 = (the space has no trace). This points to an essential difference between the properties of the classical Laplacian and the fractional Laplacian for small values of .

Subcritical boundary behaviour in average terms
We have an extension of the result for the classical Laplacian on averaged convergence to the boundary, see [36]. (a) If > − 1∕2, Proof. Assume that ⩾ 0. Let us start from (a). It is clear that, by duality, We decompose this last integral into two Using where we have used Theorem 3.4. As a consequence again by dominated convergence. The proof of (b) is analogous by using (3.5).
Let us now consider (c). As above, by duality, For the first integral, we use (3.2) with = − to deduce up to constants not depending on . For the second one we use (3.3) and (3.4) which give and therefore This completes the proof for ⩾ 0.
If changes sign, then can we apply the result have to + and − . □

Sharp weighted spaces for the Green operator
The computations above allow to complement the analysis carried out in Theorem 2.1, and improve the estimate for the optimal data from 1 loc to a weighted space. It follows the general philosophy that, due to (2.10), for any ∈ (Ω, ) we have  ∶ 1 (Ω, ( )) → 1 (Ω, ).
The result is as follows.
Proof. Take ∈ 1 (Ω, ). Then As > − − 1 by assumption, we can apply Theorem 3.4. Since > − 2 , then ( ) ≍ and, therefore, This completes the proof. □ Proof. In the notations of Theorem 3.9, note that, if < 2 , then = 0 is an admissible choice. □ For ∈ ∞ (Ω), we have shown that ( ) ≍ . To study the sharp boundary behaviour, we want to study ( )∕ . For this reason, we introduce the following definition. To prove sharp boundary behaviour we assume that the Green kernel has a -normal derivative .
• For the RFL, it follows from the Boundary Harnack Principle [5] and the boundary regularity of solutions on smooth domains [37]; for ∈ Ω, ∈ ∞ (Ω) fixed and = , we have Both factors lie in (Ω ⧵ ( )), at least for , > 0 small enough. Indeed, the first one is due to [5, Theorem 1] and is a consequence of the -harmonicity of the two involved functions close to the boundary; the second factor, instead, is more related to the smoothness of the boundary and a more classical Schauder regularity, see [37,Theorem 1.2]. The kernel has been first introduced in [1], although it is strongly related to the Martin kernel, see, for example, [6].
• For the SFL, the well-definition of is contained in [2,Lemma 14]; in this case, the proof relies on a computation on an explicit representation formula for in terms of the classical Dirichlet heat kernel.  Proof. We write Let ∈ Ω, ( ) ∈ℕ ⊂ Ω such that → as ↑ ∞, and = supp( ) ⋐ Ω. Then, up to constants, for sufficiently large Therefore, since convergence a.e. in is given by (K5), by dominated convergence Thus, The limit is, by definition, ( ). The pointwise estimate is a consequence of (3.8). □ Remark 3.16. One needs to be careful with pathological cases that do not satisfy (K5), including .
Therefore, for every ∈ Ω, Since is a bounded function in Ω and pointwise convergent, by the dominated convergence theorem as ↑ ∞ This completes the proof. □

Limit of the interior theory
A classical approach known for the usual Laplacian to recover the non-homogeneous Dirichlet boundary problem is to concentrate all mass towards the boundary. Let us recall a classical argument in the case of the usual Laplacian on the ball of radius 1. To construct a harmonic function such that = 1 on the boundary of the ball we can proceed as follows. Taking a well-chosen sequence of radial functions such that ‖ ‖ 1 (Ω) → | Ω| but such that support of is contained in the strip {1 − 1 < < 1}, one can pass to the limit ( ) by compactness. Looking at the very-weak formulation (or weak-dual one) we can check that the limit is the desired function. We now extend this argument to our general setting, and general Ω. Furthermore, ⋆ is a solution of and is given by Proof. It is clear that supp( ) = and ‖ ‖ 1 (Ω) = | Ω|. Therefore, due to Lemma 2.9, a subsequence of ( ), ( (1) ), is weakly convergent in 1 ({ ⩾ 1}) to a function 1 . A further subsequence, ( (2) ), converges in 1 ({ ⩾ 1 2 }) to a function 2 . Iterating the process, we construct sequences ( ) and functions defined on { > 1∕ }, for every ∈ ℕ. Applying Proposition 3.18 we have that for any ∈ ∞ (Ω). Therefore, For > , using = sign( − ) { ⩾1∕ } as a test function, we check that | ⩾1∕ = . We define ⋆ ( ) = ( ) for any > 1∕ ( ). Given ∈ ∞ (Ω), we ⋆ = for any > 1∕dist(supp , Ω). Therefore, for any ∈ ∞ (Ω).
If we now consider the Green representation, we get ) .
With this representation formula, we show that all convergent subsequences share a limit, and therefore the whole sequence converges. □ Remark 4.2. In Figure 4, we show a numerical simulation of the behaviour of the approximating sequence for the case of the SFL, under different values of .

Corollary 4.3. Under the assumptions and notations of Theorem 4.1, it holds
Proof. This follows by plugging (3.8) into (4.1). Indeed, Numerical solutions of L = ( 2 ) 1+ 1 < < 2 in dimension = 1. In the limit as ↑ +∞, we recover the profile of solutions to the L-harmonic problem. We implement the schemes introduced in Figure 3. where which completes the proof. □ Remark 4.4. Note that, for > − 1 2 , ⋆ has the limit rate 2 − −1 which is not accessible to solutions of the interior problem.

The L-harmonic problem
For a self-adjoint operator in our class of study it makes sense to consider the following boundary problem for some suitable test functions . In the case of the usual Laplacian, this is the non-homogeneous Dirichlet problem with data ℎ. This very weak formulation was first studied in [12]. Passing to our weak-dual formulation, (4.3) is written for any ∈ ∞ (Ω). Heuristically, (4.4) can be read as the L-harmonicity of in Ω, that is, L = 0 in Ω.
To understand this weak-dual problem, we proceed informally. If one takes = , the Dirac delta, we obtain We will see in Theorem 4.13 that Hence, in some sense (4.4) is a formulation of problem (1.8). We will devote the next paragraphs of this section to rigorously proving these intuitions.
In particular, let ⋐ Ω and take = sign( ) . Then Since this is true for all , we have that = 0 a.e. in Ω. Hence 1 = 2 . The kernel representation (4.5) follows as in the proof of Theorem 4.1, by exchanging the order of integration in (4.4). Note that this kernel representation can be rigorously justified on its own and therefore grant uniqueness. Nevertheless, since we will construct it as a limit of the interior theory, this is not needed.
We prove existence, (4.6), and (4.7) simultaneously. We split the proof into different steps. Let us first assume that 0 ⩽ ℎ ∈ ( Ω). Using the notations defined in Remark 3.1, we extend the definition of ℎ to the interior by setting ( ) = ℎ( ( )).
Let us define, for > 1 , the sequence We check that this sequence is bounded in 1 (Ω, ) by estimating We define = ( ).
We now show local 1 -weak convergence. Let ⋐ Ω. For any ⊂ we have that for some > 0, by Lemma 2.9. Therefore, the sequence is equi-integrable in and it admits a subsequence weakly convergent to some ∈ 1 ( ). That is, if we consider ∈ ∞ (Ω), with supp ⊆ , we have that On the other hand, by Proposition 3.18, we have that Therefore, for any ∈ ∞ , with supp ⊆ .
For two compacts , ′ ⋐ Ω and the corresponding , ′ built as above, we actually have = ′ in ∩ ′ . Indeed, let us consider the test function It is an admissible test function for both and ′ . Therefore, We define now We have shown above that any converging subsequence of converges weakly to over compacts. In particular, ⇀ in 1 loc . By construction solves (4.4). Passing to the limit the estimate in Theorem 2.5 we deduce that, as ↑ ∞, Moreover, in view of (2.5) we have that We deduce that the sequence converges weak-⋆ in ∞ ( ) to , and that ‖ ‖ ∞ ( ) ⩽ dist( , Ω) 2 − − ‖ℎ‖ 1 ( Ω) .
We now consider 0 ⩽ ℎ ∈ 1 ( Ω). We take an approximation sequence 0 ⩽ ℎ ∈ ( Ω) converging to ℎ in 1 ( Ω). The sequence of solutions corresponding to ℎ can be constructed through the previous step. Due to the estimates, we can pass to the limit over compacts and apply the uniqueness reasoning above to recover a function ∈ ∞ loc solution of (4.4) with data ℎ. For ℎ ∈ 1 ( Ω), we can decompose it as ℎ = ℎ + − ℎ − , construct solutions 1 and 2 corresponding to ℎ + and ℎ − and recover = 1 − 2 satisfying all properties.
as → . The above yields that
Note that ⋆ ( ) ( ) ≍ ( ) 2 −1 . (3) The case 2 − − 1 > (that is, < − 1 2 ) seems to pose problems to uniqueness. Indeed in this case It seems that problem does not have a unique solution, as ( ) + (ℎ) is also a solution for any ℎ ∈ ( Ω). Therefore, the construction of the Green operator assumed at the beginning (which chooses a single solution), seems to be made by applying some additional selection criteria. This phenomenon should be studied. (4) In trying to construct an example satisfying < − 1 2 relation, we have considered the following example: let ∈ ∞ (Ω) and consider the system where  is the Green operator of (−Δ) RFL . It seems that  is self-adjoint and, for 1∕2 < < 1, we expect its kernel to be of the form , , ∈Ω.
(5) The operators that admit exterior data = g in ℝ ⧵ Ω (for example, the RFL) have an exterior kernel, that is sometimes denoted by ℙ( , ), ∈ Ω, ∈ ℝ ⧵ Ω (in the case of the RFL it holds that ℙ( , ) = −(−Δ) ( , )). It seems reasonable that the singular solutions of type ⋆ can also be detected from the outside, as it has been done, for instance, in [6,Lemma 7] and [1, Lemma 3.6].
(6) Note that, so far, we have given all our estimates in terms of ∕ ⋆ . However, it would be nice to give an operator such that Nevertheless, the boundary behaviour of ⋆ is only known in terms of rate. An interesting question is if the following limit is defined This seems to be a further assumption on the kernel. If it is, then ≍ 1, we can set As it was pointed to us by Grubb, such operator for the RFL has been provided in [30, Corollary 6.2 and Theorem 7.1]. Also, as a matter of fact, always in the RFL case, the function defined in (5.1) is constant: this is a consequence of the equivalence between the integration by parts formula in [32,Corollary 4.5] and the one in [1, Proposition 2]. (7) In the case of the RFL, the existence of solutions which are singular at boundary can be obtained by taking the derivative of regular solutions. In Appendix A, we include an account of how positive singular solutions can be obtained, which was explained to us by Ros-Oton. It is based on a very interesting formula of chain-rule type, see (A.1). However, this argument does not seem to apply in general. In particular, it could fail in those examples where the commutation with the derivative does not hold, so that the singular rates cannot be predicted by such means. For the SFL it is easy to see that we cannot repeat the reasoning: in one dimension (say Ω = (−1, 1)), we take the first eigenfunction for the SFL ( ) = cos 2 . Then, all derivatives are bounded functions and no singularity appears. However, we have shown that the blow-up rate of the critical solution is 2( −1) . (8) It is interesting to point out that, when ≍ −2 (and it is admissible in the sense of (1.4), that is, 2 < + 1), then our main result Theorem 3.4 says that ( ) ≍ 1.
Given g ∈ ∞ ( Ω) it is therefore natural to ask whether there exists a function such that ( )( ) → g( ) as → ∈ Ω.
This would amount to studying whether the non-homogeneous Dirichlet problem has solutions for g bounded. When 2 − − 1 < 0 the solution of such problem will satisfy lim → Ω ∕ ⋆ = 0 on Ω. This would indicate that = ( ). If such exists, it will never be unique since, taking 2 ∈ ∞ (Ω),̂= ( + 2 ) = + ( 2 ) will also go to g at the boundary. Hence, it seems that, for operators L with a Green kernel satisfying 2 − − 1 < 0, and g cannot be chosen independently. In fact, for compactly supported the only possible bounded g is zero, in complete contrast to the problem for the classical Laplacian. In many cases, this inverse task of finding one or several given g turns out to be simple. If, for instance, the direct operator L is given by a singular integral, then given some bounded smooth boundary data g we can extend them to the interior of Ω as a smooth functiong, and by zero outside. Then, we can take =g and compute = L in Ω explicitly, for Ω of class 1,1 . This construction is particularly enlightening when g = 1. The natural extension to the inside is When L is the RFL, the computation yields, for ∈ Ω, The last computation is a simple although technical exercise. Note that is in 1 For the SFL, we use the kernel representation and deduce, for ∈ Ω, cf. (1.9). Curiously, in the CFL (which satisfies 2 − − 1 = 0) case, the L-harmonic problem has ⋆ ≍ 1 and we get the non-homogeneous Dirichlet problem. Thus, for = 1 one trivially has for all ∈ Ω, This shows, in particular, that = 1 is CFL-harmonic. (9) In this paper, we present optimal function data (and, under additional assumptions, measures). It is known that the RFL, SFL, and CFL not only improve integrability, but also differentiability, that is, solutions for data are in Sobolev spaces , for some > 0 and ⩾ . Then, it is natural to consider distributional data ∈ − , in the weak-dual formulation. Some results in this direction are due to Grubb [29,30]. Our techniques do not extend to distributional solutions, since we have used a theory of non-negative data: in particular, we make use of the positive part + which is not defined on distributions. This is not surprising, for our assumptions on the kernel are upper and lower bounds (and possibly continuity) but never differentiability. (10) It would interesting to see under which conditions the results that hold for ℎ ∈ 1 ( Ω) (notably, Theorems 4.6 and 4.15) can be extended to ℎ ∈ ( Ω). (11) Finally, we would like to stress that we have dealt here only with linear equations. In a nonlinear setting the situation may be even richer, as boundary singularities could be generated by the nonlinearity as remarked, for example, in [1] for the RFL and in [2] for the SFL: typically this type of behaviour is not captured by a Green-Martin representation as the one in (K0) and (4.9), because the boundary blow-up rate of these solutions exceeds the one of the reference function * as built in Theorem 4.1.

APPENDIX A: DERIVATION OF AN EXPLICIT SINGULAR SOLUTION BY X. ROS-OTON
This section is the result of conversations with Prof. X. Ros-Oton who had indicated to the authors that, at least in the RFL case, some singular solutions could be obtained by differentiating the continuous solutions of the standard theory as developed in [38], see also comment number 7 in Section 5. The objection we made that the solutions obtained by plain differentiation may change sign was taken into account. In this section, for the sake of clarity we drop the subindex RFL: (−Δ) = (−Δ) RFL . This is the way his argument proceeds in four steps. Working in arbitrary dimension, we consider functions ∈ 1,1 (ℝ ). First, we need the interesting identity This identity is proved in [38, Proof of Lemma 5.1] by calculations with integrals, but we suggest the reader to do it as an exercise, by taking Fourier transforms and manipulating the resulting formula.
Then we need the result by Getoor (see [27, section 3] or [4]) that applies to a very particular continuous solution [24,25]). The next step is new and goes as follows. If we call ( ) = (1 − | | 2 ) + and put = ⋅ ∇ we conclude that Finally, we consider = 2 − and get (−Δ) = 0 in 1 . Moreover, and this is a singular solution up to a harmless constant. This is precisely the singular solution provided by Hmissi [33] and Bogdan [6]. Using the co-area formula (see, for example, [26]), we have that ⩽ + +1 |ln |.
Let us now give a close look at (B.2). The inner ( − 1)-dimensional integral is uniformly bounded in whenever > 1 2 , then the integration in is elementary. If < 1∕2, reasoning as above we get and, in case = 1∕2, The above proves the first claim in the statement. Mind now that in the case < − 1 2 we could simply estimate • For ∈ Ω 1 we use that • For ∈ Ω 2 we use that ≍ ( ) ∫ Ω 2 ( ) + ≍ + +1 ( ) .
• For ∈ Ω 3 we use that which, rephrased, means • For the integration in Ω 4 , we use that, for ∈ Ω 4 ,