Regularity of extremal solutions of nonlocal elliptic systems

We examine regularity of the extremal solution of nonlinear nonlocal eigenvalue problems with an integro-differential operator, including the fractional Laplacian, of the form of $$ \mathcal L(u (x))= \lim_{\epsilon\to 0} \int_{\mathbb R^n\setminus B_\epsilon(x) } [u(x) - u(z)] J(z-x) dz , $$ when $J$ is a nonnegative measurable even jump kernel. We first establish stability inequalities for minimal solutions of the above system for a general nonlinearity and a general kernel. In particular, we study systems with exponential and power-type nonlinearities leading to the Gelfand system and the Lane-Emden system, with both positive and negative exponents. We consider jump kernels of the form of $J(y)=\frac{a(y/|y|)}{|y|^{n+2s}}$ where $s\in (0,1)$ and $a$ is any nonnegative even measurable function in $L^1(\mathbb {S}^{n-1})$ that satisfies ellipticity assumptions. Note that for the case of constant $a$ the above operator is called the fractional Laplacian operator and it is associated with a local equation known as the Caffaralli-Silvestre extension problem. However, in the general setting that is for all functions $a$, such an extension problem is not known. We prove regularity of extremal solutions in dimensions $n<10s$ and $ n<2s+\frac{4s}{p\mp 1}[p+\sqrt{p(p\mp1)}]$ for the Gelfand and Lane-Emden systems when $p>1$ (with positive and negative exponents), respectively. We establish the regularity results via technical integral estimates, inspired originally by Crandall-Rabinowitz, and we do not apply any local extension problem arguments. Moreover, we consider gradient systems with general nonlinearities and we establish regularity of extremal solutions in dimensions $n<4s$. As far as we know, this is the first regularity result on extremal solutions of nonlocal systems.


Introduction
Let Ω ⊂ R n be a bounded smooth domain. Consider the nonlinear nonlocal eigenvalue problem (P ) λ,γ    Lu = λF (u, v) in Ω, Lv = γG (u, v) in Ω, u, v = 0 on R n \ Ω, where λ, γ are positive parameters, F, G are smooth functions and the operator L is an integral operator of convolution type Here, J is a nonnegative measurable even jump kernel such that (1.2) R n min{|y| 2 , 1}J(y)dy < ∞.
The above nonlocal operator with a measurable kernel when the function c(x, z) is bounded between two positive constants, 0 < c 1 ≤ c 2 , is studied extensively in the literature from both theory of partial differential equations and theory of probability points of view, see the book of Bass [1] and references therein. Integro-differential equations and systems, of the above form, arise naturally in the study of stochastic processes with jumps, and more precisely in Lévy processes. A Lévy process is a stochastic process with independent and stationary increments. A special class of such processes is the so called stable processes. These are the processes that satisfy self-similarity properties, and they are 1 also the ones appearing in the Generalized Central Limit Theorem. We refer interested readers to the book of Bertoin [2] for more information. The infinitesimal generator of any isotropically symmetric stable Lévy process in R n is where µ is any nonnegative and finite measure on the unit sphere S n−1 called the spectral measure and s ∈ (0, 1). When the spectral measure is absolutely continuous, dµ(θ) = a(θ)dθ, the above operators can be rewritten in the form of (1.4) L(u(x)) = lim ǫ→0 R n \Bǫ(x) [u(x + y) + u(x − y) − 2u(x)] a(y/|y|) |y| n+2s dy, where s ∈ (0, 1) and a is any nonnegative even function in L 1 (S n−1 ). Note that the fractional Laplacian operator L = (−∆) s with 0 < s < 1 that is where c n,s is a positive constant, is the simplest stable Lévy process for dµ(θ) = c n,s dθ. Note that the above operator can be written in the form of (1.1) due to the fact that a is even. The regularity of solutions for equation Lu = f has been studied thoroughly in the literature by many experts and in this regard we refer interested to [1,9,24,36,39] and references therein. The most common assumption on the jump kernel in this context is 0 < c 1 ≤ a(θ) ≤ c 2 in S n−1 and occasionally a(θ) ≥ c 1 > 0 in a subset of S n−1 with positive measure. In this article, we consider the ellipticity assumption on the operator L of the form (1.6) 0 < c 1 ≤ inf ν∈S n−1 S n−1 |ν · θ| 2s a(θ)dθ and 0 ≤ a(θ) < c 2 for all θ ∈ S n−1 , where c 1 and c 2 are constants. Note that regularity results (interior and boundary) under such an assumption on general operator L is studied in the literature, and in this regard we refer interested readers to [36,37] and references therein. For particular nonlinearities of F and G, we consider the following Gelfand system in Ω, u, v = 0 on R n \ Ω, and the Lane-Emden system, when p > 1 in Ω, u, v = 0 on R n \ Ω, and the Lane-Emden system with singular nonlinearity, for p > 1 and when 0 < u, v < 1 Note that for the case of p = 2 the above singular nonlinearity and system is known as the MicroElec-troMechanical Systems (MEMS), see [18,28] and references therein for the mathematical analysis of such equations. In addition, we study the following gradient system with more general nonlinearities and the gradient system in Ω, u, v = 0 on R n \ Ω.
The nonlinearities f and g will satisfy various properties but will always at least satisfy (1.7) f is smooth, increasing and convex with f (0) = 1 and f superlinear at infinity.

2
A bounded weak solution pair (u, v) is called a classical solution when both components u, v are regular in the interior of Ω and (P ) λ,γ holds. Given a nonlinearity f which satisfies (1.7), the following nonlinear eigenvalue problem is now well-understood whenever s = 1 and 0 < s < 1. Brezis and Vázquez in [4] raised the question of determining the boundedness of u * for general nonlinearities f satisfying (1.7). See, for instance, [3-6, 8, 10, 14, 33-35, 41] for both local and nonlocal cases. It is known that there exists a critical parameter λ * ∈ (0, ∞), called the extremal parameter, such that for all 0 < λ < λ * there exists a smooth, minimal solution u λ of (S) λ . Here the minimal solution means in the pointwise sense. In addition for each x ∈ Ω the map λ → u λ (x) is increasing in (0, λ * ). This allows one to define the pointwise limit u * (x) := lim λրλ * u λ (x) which can be shown to be a weak solution, in a suitably defined sense, of (S) λ * . For this reason u * is called the extremal solution. It is also known that for λ > λ * , there are no weak solutions of (S) λ . The regularity of the extremal solution has been of great interests in the literature. It is known that the extremal solution can be a classical solution or it can be a singular weak solution. For the case of local second order semilinear elliptic equations that s = 1, Nedev in [34] proved that u * is bounded if general nonlinearity f satisfies (1.7) and n ≤ 3. This was extended to fourth dimensions when Ω is a convex domain by Cabré in [5]. For the particular nonlinearity f (u) = e u , known as the Gelfand equation, the regularity is shown u * ∈ L ∞ (Ω) for dimensions n < 10 by Crandall and Rabinowitz in [14], see also [20]. If Ω is a radial domain in R n with n < 10 the regularity is shown in [6] when f is a general nonlinearity satisfying conditions (1.7) but without the convexity assumption. In view of the above result for the exponential nonlinearity, this is optimal. Note that for the case of Ω = B 1 , the classification of all radial solutions to this problem was originally done by Liouville in [30] for n = 2 and then in higher dimensions in [27, 31,33] and references therein. For power nonlinearity f (u) = (1 + u) p for p > 1, known as the Lane-Emden equation, the regularity is shown for dimensions by Farina in [19]. The critical dimension and critical exponent can be rephrased as 1 < p < p c (n) when p c is called Joseph-Lundgren [27] exponent and given by if n ≥ 11, (1.8) Note that Joseph-Lundgren exponent is larger than the Sobolev exponent. In addition for the case of singular nonlinearity f (u) = 1 (1−u) 2 , known as the MEMS equation, the regularity is proved for n ≤ 7 by Ghoussoub et al. in [17,18,25]. Note the latter is extended to the case of f (u) = 1 (1−u) p for p > 1 for dimensions Some of the above results are given for the case of fractional Laplacian 0 < s < 1 by Ros-Oton and Serra in [35] and also [38]. Note also that it is well-known that there is a correspondence between the regularity of stable solutions on bounded domains and the Liouville theorems for stable solutions on the entire space, via rescaling and a blow-up procedure. In this regard classifications of solutions of the above equation on the entire space is studied in the literature and we refer interested readers to [11, 15, 18-20, 23, 29] and references therein. For the case of systems, as discussed in [13,21,32], set Q := {(λ, γ) : λ, γ > 0} and define (1.9) U := {(λ, γ) ∈ Q : there exists a smooth solution (u, v) of (P ) λ,γ } .
We assume that F (0, 0), G(0, 0) > 0. A simple argument shows that if F is superlinear at infinity for u, uniformly in v, then the set of λ in U is bounded. Similarly we assume that G is superlinear at infinity for v, uniformly in u, and hence we get U is bounded. We also assume that F, G are increasing in each variable. This allows the use of a sub/supersolution approach and one easily sees that if (λ, γ) ∈ U then so is (0, λ] × (0, γ]. One also sees that U is nonempty. We now define Υ := ∂U ∩ Q, which plays the role of the extremal parameter λ * . Various properties of Υ are known, see [32]. Given (λ * , γ * ) ∈ Υ set σ := γ * λ * ∈ (0, ∞) and define We let (u λ , v λ ) denote the minimal solution (P ) λ,σλ for λ * 2 < λ < λ * . One easily sees that for each x ∈ Ω that u λ (x), v λ (x) are increasing in λ and hence we define and we call (u * , v * ) the extremal solution associated with (λ * , γ * ) ∈ Υ. Under some very minor growth assumptions on F and G one can show that (u * , v * ) is a weak solution of (P ) λ * ,γ * . For the rest of this article we refer to (u * , v * ) as (u, v). For the case of local Laplacian operator, Cowan and the author in [13] proved that the extremal solution of (H) λ,γ when Ω is a convex domain is regular provided 1 ≤ n ≤ 3 for general nonlinearities f, g ∈ C 1 (R) that satisfying (1.7). This can be seen as a counterpart of the Nedev's result for elliptic gradient systems. For radial solutions, it is also shown in [13] that stable solutions are regular in dimensions 1 ≤ n < 10 for general nonlinearities. This is a counterpart of the regularity result of Cabré-Cappella [6] and Villegas [41] for elliptic gradient systems. For the local Gelfand system, regularity of the extremal solutions is given by Cowan in [12] and by Dupaigne et al. in [16] when n < 10.
In the current article, we provide regularity of the extremal solutions of nonlocal Gelfand, Lane-Emden and MEMS systems for the operator L given by (1.4). In addition, we provide a counterpart of the Nedev regularity result for system (H) λ,γ . To do so, we first establish integral estimates for minimal solutions of systems with a general nonlocal operator L given by (1.1) when J is a nonnegative measurable even jump kernel. Then, we apply nonlocal Sobolev embedding arguments to conclude boundedness of the extremal solutions. Note that when λ = γ the systems (G) λ,γ , (E) λ,γ and (M ) λ,γ with particular nonlinearities turn into scalar equations. As it was stated earlier, there is a correspondence between the regularity of stable solutions on bounded domains and the Liouville theorems for stable solutions on the entire space, via rescaling and a blow-up procedure. This enables one to establish Liouville theorems for the corresponding systems.
Here is how this article is structured. In Section 2, we provide regularity theory for nonlocal operators. In Section 3, we establish various stability inequalities for minimal solutions of systems with a general nonlocal operator of the form (1.1) with a nonnegative measurable even jump kernel. In Section 4, we provide some technical integral estimates for stable solutions of systems introduced in the above. In Section 5, we apply the integral estimates to establish regularity of extremal solutions for nonlocal Gelfand, Lane-Emden and MEMS systems with exponential and power nonlinearities. In addition, we provide regularity of extremal solutions for the gradient system (H) λ,γ with general nonlinearities and also for particular power nonlinearities.

Preliminaries
In this section, we provide regularity results not only to the fractional Laplacian, but also to more general integro-differential equations. We omit the proofs in this section and refer interested readers to corresponding references. Let us start with the following classical regularity result concerning embeddings for the Riesz potential, see the book of Stein [40].
Theorem 2.1. Suppose that 0 < s < 1, n > 2s and f and u satisfy in the sense that u is the Riesz potential of order 2s of f . Let u, f ∈ L p (R n ) when 1 ≤ p < ∞.
(i) For p = 1, there exists a positive constant C such that (ii) For 1 < p < n 2s , there exists a positive constant C such that where [·] C β (R n ) denotes the C β seminorm. Here the constant C depending only on n, s and p.
The above theorem is applied by Ros-Oton and Serra in [35] to establish the following regularity theory for the fractional Laplacian. See also [36] for the boundary regularity results.
(iii) For n 2s < p < ∞, there exists a positive constant C such that Here the constant C depending only on n, s, p, r and Ω.
In what follows we provide a counterpart of the above regularity result for general integro-differential operators given by (1.4). These operators are infinitessimal generators of stable and symmetric Lévy processes and they are uniquely determined by a finite measure on the unit sphere S n−1 , often referred as the spectral measure of the process. When this measure is absolutely continuous, symmetric stable processes have generators of the form (1.4) where 0 < s < 1 and a is any nonnegative function L 1 (S n−1 ) satisfying a(θ) = a(−θ) for θ ∈ S n−1 . The regularity theory for general operators of the form (1.4) has been recently developed by Fernández-Real and Ros-Oton in [24]. In order to prove this result, authors apply results of [26] to study the fundamental solution associated to the operator L in view of the one of the fractional Laplacian.

Proposition 2.2.
Let Ω ⊂ R n be any bounded domain, 0 < s < 1 and f ∈ L 2 (Ω). Let u be any weak solution of where the operator L is given by (1.4) and the ellipticity condition (1.6) holds. Assume that f ∈ L r (Ω) for some r.
(iii) For n 2s < r < ∞, there exists a positive constant C such that Here the constant C depending only on n, s, r, Ω and ellipticity constants.
We end this section with this point that (u, v) is a weak solution of (P ) λ, when ζ, η and ζ, η are bounded in Ω and ζ, η ≡ 0 on ∂Ω. Any bounded weak solution is a classical solution, in the sense that it is regular in the interior of Ω, continuous up to the boundary, and (P ) λ,γ holds pointwise. Note that for the case of local operators, that is s = 1, the above notion of weak solution is consistent with the one introduced by Brezis et al. in [3,4].

Stability Inequality
In this section, we provide stability inequalities for minimal solutions of system (P ) λ,γ for various nonlinearities F and G. We start with the following technical lemma in regards to nonlocal operator L with even symmetric kernel J.
Lemma 3.1. Assume that an operator L is given by (1.5) with a measurable symmetric kernel J(x, z) = J(x − z) that is even. Then, where f, g ∈ C 1 (R n ) and the integrals are finite.
Proof. The proof is elementary and we omit it here.
We now establish a stability inequality for minimal solutions of system (P ) λ,γ . Note that for the case of local operators this inequality is established by the author and Cowan in [13] and in [21].
Assume that J is a measurable even kernel and F v G u ≥ 0. Then, for test functions ζ, η so that ζ, η = 0 in R n \ Ω.
Proof. Since u λ , v λ are increasing in λ, differentiating (P ) λ,γ with respect to λ we get Note that the following lower-bound holds for the left-hand side of the above equality Integrating the above we end up with Applying Lemma 3.1, we have Note that for a, b, c, d ∈ R when ab < 0 we have Since each ∂ λ u λ does not change sign, we have in the above inequality and from the fact that ab = −∂ λ u λ (x)∂ λ u λ (z) < 0, we conclude Therefore, This together with (3.2) complete the proof.
Following ideas provided in the above, we provide stability inequalities for minimal solutions of Gelfand, Lane-Emden and MEMS systems with exponential and power-type nonlinearities.
for test functions ζ so that ζ = 0 in R n \ Ω.

Integral Estimates for Stable Solutions
In this section, we establish some technical integral estimates for stable solutions of systems. Most of the ideas and methods applied in this section are inspired by the ones developed in the literature, see for example [16,19,20,22,23]. We start with the Gelfand system. Proof. Multiply the second equation of (G) λ,γ with e u − 1 and integrate to get Note that for a, b ∈ R, one can see that Applying the above inequality for a = u(x) and b = u(z), we obtain From the above, we conclude Test the stability inequality, Corollary 3.1, on ζ = e u 2 − 1 to get Combining above inequalities, we conclude Applying the Young's inequality e u/2 ≤ e u 4 + 1, we conclude we complete the proof.
We now provide a counterpart of the above estimate for stable solutions of (E) λ,γ .
is a solution of (E) λ,γ when the associated stability inequality (3.4) holds.
Then, there exists a positive constant C λ,γ,|Ω| = C(λ, γ, |Ω|) such that Proof. Multiply the second equation of (E) λ,γ with (1 + u) p − 1 and integrate to get From Lemma 3.1, we get Note that for a, b ∈ R, one can see that Applying the above inequality for a = u(x) and b = u(z), we obtain From the above, we conclude Test the stability inequality, Corollary 3.1, on ζ = (1 + u) Combining above inequalities, we conclude Expanding the left-hand side of the inequality and rearranging we get where we have used the inequality a ≤ ǫ 2 a 2 + 1 2ǫ for any ǫ > 0. Similarly, Note that from the Cauchy-Schwarz inequality we get From this and multiplying both sides of the above (4.8) and (4.9) we conclude for small ǫ > 0. Note that p 2 − (p+1) 2 4p 2 > 0 when p > 1. Therefore, taking small enough ǫ > 0 and applying the Hölder's inequality, we complete the proof.
Here is a counterpart of the above estimate for stable solutions of (M ) λ,γ .   Proof. The proof is similar to the one provide in Lemma 4.2. Multiply the second equation of (M ) λ,γ with (1 − u) −p − 1 and integrate to get Note that for a, b ∈ R, one can see that Applying the above inequality for a = u(x) and b = u(z), we obtain From the above, we conclude Test the stability inequality, Corollary 3.1, on ζ Combining above inequalities, we conclude Expanding the left-hand side of the inequality and rearranging we get where we have used the inequality a ≤ ǫ 2 a 2 + 1 2ǫ for any ǫ > 0. Similarly, Note that from the Cauchy-Schwarz inequality we get From this and multiplying both sides of the above (4.8) and (4.9) we conclude for small ǫ > 0. Note that p 2 − ( (p−1) 2 4p ) 2 > 0 when p > 0. Therefore, taking small enough ǫ > 0 and applying the Hölder's inequality when p > 1, we complete the proof.
In the next lemmata, we provide integral L q (Ω) estimates for Gelfand, Lane-Emden and MEMS systems. We start with the Gelfand system and establish a relation between Ω e where C ǫ,λ,γ,|Ω| is a positive constant.
Proof. Multiply the second equation of (G) λ,γ with e tu − 1 when t > 1 2 is a constant. Integrating implies that λ Note that for a, b ∈ R, one can see that Applying the above inequality for a = tu(x) and b = tu(z), we obtain From the above, we conclude Test the stability inequality, Corollary 3.1, on ζ = e tu 2 − 1 to get Combining above inequalities, we conclude (4.11) λγ On the other hand, from the Young inequality we have (4.12) where ǫ is a positive constant. In addition, from the Hölder inequality we get Combining (4.12), (4.13) and (4.14) proves the first inequality in (4.10). With similar arguments one can show the second inequality.
We now consider the Lane-Emden system and establish a relation between Ω (1 + u) where C ǫ,λ,γ,|Ω| is a positive constant.
Proof. Let t > 1 be a constant. Multiply the second equation of (E) λ,γ with (1 + u) t − 1 and integrate to get Note that for a, b ∈ R, one can see that Applying the above inequality for a = u(x) and b = u(z), we obtain From the above, we conclude Test the stability inequality, Corollary 3.1, on ζ = (1 + u) Combining above inequalities, we conclude Expanding the left-hand side of the inequality and rearranging we get where we have used the inequality a ≤ ǫ 2 a 2 + 1 2ǫ for any ǫ > 0. From the Hölder's inequality we get where β = 2(t+1) 2t−p+1 . This and (4.20) completes the proof of the first estimate in (4.15). Similarly, one can show the second estimate.
We now consider the MEMS system with singular power nonlinearities and establish a relation between for t > 1. Then, for some constant 0 < ǫ < 1 we get where C ǫ,λ,γ,|Ω| is a positive constant.
Proof. Let t > 1 be a constant. Multiply the second equation of (M ) λ,γ with (1 − u) −t − 1 and integrate to get Note that for a, b ∈ R, one can see that Applying the above inequality for a = u(x) and b = u(z), we obtain From the above, we conclude Test the stability inequality, Corollary 3.1, Combining above inequalities, we conclude (4.19) λγp 13 Expanding the left-hand side of the inequality and rearranging we get where we have used the inequality a ≤ ǫ 2 a 2 + 1 2ǫ for any ǫ > 0 and a ∈ R. From the Hölder's inequality we get This and (4.20) completes the proof of the first estimate in (4.15). Similarly, one can show the second estimate.
In regards to the gradient system with superlinear nonlinearities satisfying (1.7) we establish an integral estimate that yields L 2 (Ω) of the function f ′ (u)g ′ (v). We then use this to conclude estimates on the nonlinearities of the gradient system. Our methods and ideas in the proof are inspired by the ones developed in [13] and originally by Nedev in [34].  Let (λ * , γ * ) ∈ Υ and (u, v) denote the extremal solution associated with (H) λ * ,γ * . Then, there exists a positive constant C < ∞ such that Proof. We obtain uniform estimates for any minimal solution (u, v) of (H) λ,γ on the ray Γ σ and then one sends λ ր λ * to obtain the same estimate for (u * , v * ). Let (u, v) denote a smooth minimal solution of (H) λ,γ on the ray Γ σ and put ζ := f ′ (u) − a and η : Note that for all a, b ∈ R, one can see that when h 1 (s) := s 0 |f ′′ (w)| 2 dw. Similar inequality holds for the function g that is 14 when h 2 (s) := s 0 |g ′′ (w)| 2 dw. Set a = u(x) and b = u(z) and a = v(x) and b = v(z) in the above inequalities to conclude From the equation of system and Lemma 3.1, we get From this and we conclude that Given the assumptions, there is some M > 1 large enough and 0 < δ < 1 that for all u ≥ M we have We now estimate the last integral in the above. Let k ≥ 1 denote a natural number. Then, We now discuss the last integral in the above. Note that this integral is bounded above by From the above estimates, we conclude that for sufficiently large M and for all 1 ≤ k there is some positive constant C(k, M ) and 0 < δ < 1 that Applying the same argument, one can show that for sufficiently large M and for all 1 ≤ k there is some positive constant D(k, M ) and 0 < ǫ < 1 that Note that f ′′ (u), g ′′ (v) → ∞ when u, v → ∞. This implies that Now set k to be sufficiently large and substitute (4.23) and (4.24) in (4.22) to complete the proof. Note that see that all the integrals in (4.22) are bounded independent of λ and γ.

Regularity of the Extremal Solutions
In this section, we apply the integral estimates established in the latter section to prove regularity results for extremal solutions of systems mentioned in the introduction earlier. We shall start with the systems with particular nonlinearities. Note that for the case of Laplacian operator L = −∆, that is s = 1, the following theorem recovers the optimal dimensions provided for associated local systems.
Even though the above theorem is optimal as s → 1, it is not optimal for smaller values of 0 < s < 1. In this regard, consider the case of λ = γ and the Gelfand system turns into (−∆) s u = λe u in the entire space .
In particular, the extremal solution should be bounded when n ≤ 7 for 0 < s < 1. We refer interested to [35,38] for more details. Now, consider the case of Lane-Emden equation (−∆) s u = λu p in the entire space R n . It is also known that the explicit singular solution u s (x) = A|x| − 2s p−1 where the constant A is given by .
As s → 1, the above inequality is consistent with the dimensions given in (5.2). For more information, we refer interested readers to Wei and the author in [23] and Davila et al. in [15]. Given above, proof of the optimal dimension for regularity of extremal solutions remains an open problem. We now establish a regularity result for the gradient system (H) λ,γ with general nonlinearities f and g satisfying certain assumptions. Note that for the case of local scalar equations such results are provided by Nedev in [34] for n = 3 and Cabré in [5] for n = 4. For the case of scalar equation with the fractional Laplacian operator, Ros-Oton and Serra in [35] established regularity results for dimensions n < 4s when 0 < s < 1. For the case of local gradient systems, that is when s = 1, such a regularity result is established by the author and Cowan in [13] in dimensions n ≤ 3.