Eigenfunctions and Fundamental Solutions of the Fractional Two-Parameter Laplacian

Abstract

We deal with the following fractional generalization of the Laplace equation for rectangular domains (𝑥,𝑦)∈(𝑥₀,𝑋₀)×(𝑦₀,𝑌₀)⊂ℝ₊×ℝ₊, which is associated with the Riemann-Liouville fractional derivatives Δ^𝛼,𝛽𝑢(𝑥,𝑦)=𝜆𝑢(𝑥,𝑦), Δ^𝛼,𝛽∶=𝐷_𝑥^1+𝛼₀₊+𝐷_𝑦^1+𝛽₀₊, where 𝜆∈ℂ, (𝛼,𝛽)∈[0,1]×[0,1]. Reducing the left-hand side of this equation to the sum of fractional integrals by 𝑥 and 𝑦, we then use the operational technique for the conventional right-sided Laplace transformation and its extension to generalized functions to describe a complete family of eigenfunctions and fundamental solutions of the operator Δ^𝛼,𝛽 in classes of functions represented by the left-sided fractional integral of a summable function or just admitting a summable fractional derivative. A symbolic operational form of the solutions in terms of the Mittag-Leffler functions is exhibited. The case of the separation of variables is also considered. An analog of the fractional logarithmic solution is presented. Classical particular cases of solutions are demonstrated.

1. Introduction

Let 𝐷^𝛾_𝑎+𝑓 and 𝐼^𝛾_𝑎+𝑓 be the Riemann-Liouville fractional derivative and integral of order 𝛾>0 defined by [1, 2]

(𝐷^𝛾_𝑎+𝑑𝑓)(𝑥)=()𝑑𝑥^𝑛1Γ(𝑛−𝛾)^𝑥_𝑎𝑓(𝑡)(𝑥−𝑡)^{𝛾−𝑛+1}𝑑𝑡,𝑎,𝑥>0,𝑛=[𝛾]+1,(𝐼^𝛾_𝑎+1𝑓)(𝑥)=Γ(𝛾)^𝑥_𝑎𝑓(𝑡)(𝑥−𝑡)^1−𝛾𝑑𝑡,𝑎,𝑥>0, (1.1)(1.2) where [] means the integer part of 𝛾. Consider a class of the linear nonhomogeneous differential equations:

(𝐷_𝑥^1+𝛼₀₊𝑢)(𝑥,𝑦)+(𝐷_𝑦^1+𝛽₀₊𝑢)(𝑥,𝑦)−𝜆𝑢(𝑥,𝑦)=𝑓(𝑥,𝑦), (1.3) where 𝜆∈ℂ, (𝛼,𝛽)∈[0,1]×[0,1], 𝑓(𝑥,𝑦), 0≤𝑥₀≤𝑥≤𝑋₀<∞, 0≤𝑦₀≤𝑦≤𝑌₀<∞ is a prescribed function, and 𝑢(𝑥,𝑦) is to be determined. Denoting by

Δ^𝛼,𝛽∶=𝐷_𝑥^1+𝛼₀₊+𝐷_𝑦^1+𝛽₀₊, (1.4) equation (1.3) can be written in the form (Δ^𝛼,𝛽−𝜆𝐸)𝑢=𝑓, where 𝐸 is the identity operator. When 𝛼=𝛽=1, we come out with the classical Poisson equation. Therefore we call fractional partial differential equation (1.3) the fractional two-parameter Poisson equation (FPE). Its homogeneous analog is naturally called the fractional Laplace equation (FLE) or fractional two-parameter Laplacian.

In this paper we present a general operational approach [3] to describe eigenfunctions and fundamental solutions of the fractional two-parameter Laplacian based on the conventional right-sided Laplace transform [4]

𝐹(𝑠)=^∞_𝑇𝑓𝑓(𝑡)𝑒^−𝑠𝑡𝑑𝑡,𝑇_𝑓≥0,Re𝑠>𝑎₀ (1.5) of absolutely integrable functions 𝑓∈𝐿₁((𝑇_𝑓,∞);𝑒^−𝑎0𝑡𝑑𝑡) with respect to the measure 𝑒^−𝑎0𝑡𝑑𝑡 and its distributional analog

𝐹(𝑠)=⟨𝑓(𝑡),𝑒^−𝑠𝑡⟩ (1.6) in Zemanian's space ℒ^′(𝑎₀) defined below. Operational solutions will be written in terms of the generalized Mittag-Leffler function 𝐸_𝜇,𝜈(𝑧), [1, 2, 5] which is defined in terms of the power series:

𝐸_𝜇,𝜈(𝑧)=^∞_𝑛=0𝑧^𝑛Γ(𝜇𝑛+𝜈),𝜇>0,𝜈∈ℝ,𝑧∈ℂ. (1.7) In particular, the function 𝐸_𝜇,𝜈(𝑧) is entire of the order 𝜌=1/𝜇 and type 𝜎=1. The exponential function and trigonometric and hyperbolic functions are expressed through (1.7) as follows:

𝐸_1,1(𝑧)=𝑒^𝑧,𝐸_2,1(−𝑧²)=cos𝑧,𝐸_2,1(𝑧)=cosh𝑧,𝑧𝐸_2,2(−𝑧²)=sin𝑧,𝑧𝐸_2,2(𝑧²)=sinh𝑧. (1.8) We will consider in the sequel the existence and uniqueness of general solutions of the fractional Laplacian and its particular cases. Possible applications and an investigation of the fractional two-parameter Poisson equation (1.3) are still out of the framework of this paper and will be done in forthcoming articles of the author.

2. Eigenfunctions and Fundamental Solutions of the Fractional Laplace Equation

We begin with the following.

Definition 2.1 (see [1]). By 𝐴𝐶^𝑛([𝑎,𝑏]), 𝑛∈ℕ, one denotes the class of functions 𝑓(𝑥), which are continuously differentiable on the segment [𝑎,𝑏] up to the order 𝑛−1 and 𝑓^(𝑛−1)(𝑥) is absolutely continuous on [𝑎,𝑏].

It is known [1] that the class 𝐴C^𝑛([𝑎,𝑏]) contains only functions represented in the form

1𝑓(𝑥)=(𝑛−1)!^𝑥_𝑎(𝑥−𝑡)^𝑛−1𝜑(𝑡)𝑑𝑡+^𝑛−1_𝑘=0𝑐_𝑘(𝑥−𝑎)^𝑘, (2.1) where 𝜑(𝑡)∈𝐿₁(𝑎,𝑏) and 𝑐_𝑘 are arbitrary constants. It is not difficult to find that 𝜑(𝑡)=𝑓^(𝑛)(𝑡), 𝑐_𝑘=𝑓^(𝑘)(𝑎)/𝑘!, Moreover, if 𝑓(𝑥)∈𝐴𝐶^𝑛([𝑎,𝑏]), then fractional derivative (1.1) exists almost everywhere and can be represented by the formula

(𝐷^𝛾_𝑎+𝑓)(𝑥)=^𝑛−1_𝑘=0𝑓^(𝑘)(𝑎)Γ(1+𝑘−𝛾)(𝑥−𝑎)^𝑘−𝛾+1Γ(𝑛−𝛾)^𝑥_𝑎𝑓^(𝑛)(𝑡)(𝑥−𝑡)^{𝛾−𝑛+1}𝑑𝑡,𝑛=[𝛾]+1. (2.2)

Definition 2.2 (see [1]). By 𝐼^𝛾_𝑎+(𝐿₁) denotes the class of functions 𝑓 represented by the left-sided fractional integral (1.2) of a summable function, that is, 𝑓=𝐼^𝛾_𝑎+𝜑, 𝜑∈𝐿₁(𝑎,𝑏).

A description of this class is given by the following.

Theorem 2.3 (see [1]). A function 𝑓(𝑥)∈𝐼^𝛾_𝑎+(𝐿₁), 𝛾>0 if and only if (𝐼^𝑛−𝛾_𝑎+𝑓)(𝑥)∈𝐴𝐶^𝑛([𝑎,𝑏]), 𝑛=[𝛾]+1 and (𝐼^𝑛−𝛾_𝑎+𝑓)^(𝑘)(𝑎)=0, 𝑘=0,1,…,𝑛−1.

Definition 2.4 (see [1]). One will say that a function 𝑓∈𝐿₁(𝑎,𝑏) has a summable fractional derivative (𝐷^𝛾_𝑎+𝑓)(𝑥) if (𝐼^𝑛−𝛾_𝑎+𝑓)(𝑥)∈𝐴𝐶^𝑛([𝑎,𝑏]), 𝑛=[𝛾]+1.

If (𝐷^𝛾_𝑎+𝑓)(𝑥)=(𝑑/𝑑𝑥)^𝑛(𝐼^𝑛−𝛾_𝑎+𝑓)(𝑥) exists in the ordinary sense, that is, (𝐼^𝑛−𝛾_𝑎+𝑓)(𝑥) is differentiable in each point up to the order 𝑛, then 𝑓(𝑥) evidently admits the derivative (𝐷^𝛾_𝑎+𝑓)(𝑥) in the sense of Definition 2.4.

So, if 𝑓(𝑥)∈𝐼^𝛾_𝑎+(𝐿₁), then (𝐼^𝛾_𝑎+𝐷^𝛾_𝑎+𝑓)(𝑥)=𝑓(𝑥). Otherwise if 𝑓 just admits a summable fractional derivative, then the composition of fractional operators (1.1) and (1.2) can be written in the form (see [1])

(𝐼^𝛾_𝑎+𝐷^𝛾_𝑎+𝑓)(𝑥)=𝑓(𝑥)−^𝑛−1_𝑘=0(𝑥−𝑎)^{𝛾−𝑘−1}Γ(𝛾−𝑘)(𝐼^𝑛−𝛾_𝑎+𝑓)^{(𝑛−𝑘−1)}(𝑎),𝑛=[𝛾]+1. (2.3) Nevertheless we note that (𝐷^𝛾_𝑎+𝐼^𝛾_𝑎+𝑓)(𝑥)=𝑓(𝑥) for any summable function 𝑓.

Consider now the eigenfunction problem for the fractional Laplace equation in the rectangular domain (𝑥,𝑦)∈(𝑥₀,𝑋₀)×(𝑦₀,𝑌₀)

(𝐷_𝑥^1+𝛼₀₊𝑢)(𝑥,𝑦)+(𝐷_𝑦^1+𝛽₀₊𝑢)(𝑥,𝑦)=𝜆𝑢(𝑥,𝑦), (2.4) where 𝜆∈ℂ, (𝛼,𝛽)∈(0,1]×(0,1] in the following three cases:

(i)𝑢(𝑥,𝑦) belongs to classes 𝐼_𝑥^1+𝛼₀₊(𝐿₁), 𝐼_𝑦^1+𝛽₀₊(𝐿₁) by 𝑥∈[𝑥₀,𝑋₀] and 𝑦∈[𝑦₀,𝑌₀], respectively;(ii)𝑢(𝑥,𝑦) admits a summable fractional derivative (𝐷_𝑥^1+𝛼₀₊𝑢)(𝑥,𝑦) by 𝑥∈[𝑥₀,𝑋₀] and belongs to 𝐼_𝑦^1+𝛽₀₊(𝐿₁) by 𝑦∈[𝑦₀,𝑌₀] or vice versa;(iii)𝑢(𝑥,𝑦) admits summable fractional derivative (𝐷_𝑥^1+𝛼₀₊𝑢)(𝑥,𝑦), (𝐷_𝑦^1+𝛽₀₊𝑢)(𝑥,𝑦) by 𝑥∈[𝑥₀,𝑋₀] and 𝑦∈[𝑦₀,𝑌₀], respectively.

Theorem 2.5. In case (i) trivial solution of (2.4) is the only solution.

Proof. Indeed, taking the operator 𝐼_𝑥^1+𝛼₀₊ from both sides of (2.4) and using the identity (𝐼_𝑥^1+𝛼₀₊𝐷_𝑥^1+𝛼₀₊𝑢)(𝑥,𝑦)=𝑢(𝑥,𝑦) it becomes 𝑢(𝑥,𝑦)+(𝐼_𝑥^1+𝛼₀₊𝐷_𝑦^1+𝛽₀₊𝑢)(𝑥,𝑦)−𝜆(𝐼_𝑥^1+𝛼₀₊𝑢)(𝑥,𝑦)=0. (2.5) Hence, applying the operator 𝐼_𝑦^1+𝛽₀₊ to both sides of (2.5), we use the fact that due Fubini's theorem this operator commutes with 𝐼_𝑥^1+𝛼₀₊. Then we obtain (𝐼_𝑥^1+𝛼₀₊𝑢)(𝑥,𝑦)+(𝐼_𝑦^1+𝛽₀₊𝑢)(𝑥,𝑦)−𝜆(𝐼_𝑥^1+𝛼₀₊𝐼_𝑦^1+𝛽₀₊𝑢)(𝑥,𝑦)=0. (2.6) Hence from conditions of the theorem we observe that fractional integrals of the equation (2.6) are Laplace-transformable functions. Therefore we may act on (2.6) by the conventional right-sided Laplace transform (1.5), let say, by 𝑦 with 𝑇_𝑓=𝑦₀. Taking into account its convolution and operational properties [3] after straightforward calculations we arrive at the following second kind homogeneous integral equation of the Volterra type: 𝑈(𝑥,𝑠)+(𝑠^1+𝛽−𝜆)Γ(1+𝛼)^𝑥_𝑥0(𝑥−𝑡)^𝛼𝑈(𝑡,𝑠)𝑑𝑡=0, (2.7) where 𝜆∈ℂ, 𝑠∈ℂ, Re𝑠≥𝑎₀>0, 𝑠^1+𝛽=|𝑠|^1+𝛽𝑒^𝑖𝜃,𝜃=arg𝑠∈(−𝜋/2,𝜋/2) and 𝑈(𝑥,𝑠)=^𝑌0_𝑦0𝑒^−𝑠𝑡𝑢(𝑥,𝑡)𝑑𝑡. (2.8) Appealing to [5, Chapter 3] we find that (2.7) has the only trivial solution in the space of summable functions and 𝑈(𝑥,𝑠)∈𝐿₁(𝑥₀,𝑋₀) because 𝑈(𝑥,𝑠)∈𝐼_𝑥^1+𝛼₀₊(𝐿₁) for each 𝑠. Cancelling the Laplace transform and using its uniqueness property for summable functions we get 𝑢(𝑥,𝑦)=0. Theorem 2.5 is proved.
In case (ii), (2.6) should be substituted by the following equalities (see (2.3)):
(𝐼_𝑥^1+𝛼₀₊𝑢)(𝑥,𝑦)+(𝐼_𝑦^1+𝛽₀₊𝑢)(𝑥,𝑦)−𝜆(𝐼_𝑥^1+𝛼₀₊I_𝑦^1+𝛽₀₊=𝑢)(𝑥,𝑦)(𝑥−𝑥₀)^𝛼−1Γ(𝛼)(𝐼_𝑦^1+𝛽₀₊𝑓₀)(𝑦)+(𝑥−𝑥₀)^𝛼Γ(1+𝛼)(𝐼_𝑦^1+𝛽₀₊𝑓₁)(𝑦), (2.9) or (𝐼_𝑥^1+𝛼₀₊𝑢)(𝑥,𝑦)+(𝐼_𝑦^1+𝛽₀₊𝑢)(𝑥,𝑦)−𝜆(𝐼_𝑥^1+𝛼₀₊𝐼_𝑦^1+𝛽₀₊=𝑢)(𝑥,𝑦)(𝑦−𝑦₀)^𝛽−1Γ(𝛽)(𝐼_𝑥^1+𝛼₀₊ℎ₀)(𝑥)+(𝑦−𝑦₀)^𝛽Γ(1+𝛽)(𝐼_𝑥^1+𝛼₀₊ℎ₁)(𝑥), (2.10) where we denoted by𝑓₀(𝑦)=(𝐼_𝑥^1−𝛼₀₊𝑢)(𝑥₀𝑓,𝑦),₁(𝑦)=(𝐷^𝛼_𝑥0+𝑢)(𝑥₀ℎ,𝑦),₀(𝑥)=(𝐼_𝑦^1−𝛽₀₊𝑢)(𝑥,𝑦₀ℎ),₁(𝑥)=(𝐷^𝛽_𝑦0+𝑢)(𝑥,𝑦₀) (2.11)(2.12)(2.13)(2.14) Cauchy's fractional initial conditions. Treating, for instance, (2.9) we take the Laplace transform from its both sides and arrive at the following integral equation: 𝑈(𝑥,𝑠)+(𝑠^1+𝛽−𝜆)Γ(1+𝛼)^𝑥_𝑥0(𝑥−𝑡)^𝛼𝑈(𝑡,𝑠)𝑑𝑡=𝐹(𝑥,𝑠),𝑥∈(𝑥₀,𝑋₀),𝜆∈ℂ, (2.15) where 𝐹(𝑥,𝑠)=(𝑥−𝑥₀)^𝛼−1𝐹Γ(𝛼)⁰(𝑠)+(𝑥−𝑥₀)^𝛼𝐹Γ(1+𝛼)¹𝐹(𝑠),_𝑖(𝑠)=^𝑌0_𝑦0𝑒^−𝑠𝑡𝑓_𝑖(𝑡)𝑑𝑡,𝑖=0,1. (2.16)(2.17) It is known [5] that a unique solution of (2.15) in the class of summable functions is 𝑈(𝑥,𝑠)=𝐹(𝑥,𝑠)−(𝑠^1+𝛽−𝜆)^𝑥_𝑥0(𝑥−𝑡)^𝛼𝐸_{1+𝛼,1+𝛼}(−(𝑠^1+𝛽−𝜆)(𝑥−𝑡)^1+𝛼)𝐹(𝑡,𝑠)𝑑𝑡, (2.18) which involves as the kernel the generalized Mittag-Leffler function (1.7). Next, substituting (2.16) into (2.15), using (1.7), (2.17), index laws for fractional operators [1], and the estimates |(𝑠^1+𝛽−𝜆)^𝑥_𝑥0(𝑥−𝑡)^𝛼𝐸_{1+𝛼,1+𝛼}(−(𝑠^1+𝛽−𝜆)(𝑥−𝑡)^1+𝛼)𝐹(𝑡,𝑠)𝑑𝑡|≤|𝐹₀(𝑠)|^∞_𝑛=0(|𝑠|^1+𝛽+|𝜆|)^𝑛+1Γ((𝑛+1)(1+𝛼))Γ(𝛼)^𝑥_𝑥0(𝑥−𝑡)^{𝑛(1+𝛼)+𝛼}(𝑡−𝑥₀)^𝛼−1𝑑𝑡+|𝐹₁(𝑠)|^∞_𝑛=0(|𝑠|^1+𝛽+|𝜆|)^𝑛+1Γ((𝑛+1)(1+𝛼))Γ(1+𝛼)^𝑥_𝑥0(𝑥−𝑡)^{𝑛(1+𝛼)+𝛼}(𝑡−𝑥₀)^𝛼𝑑𝑡≤(|𝑠|^1+𝛽+|𝜆|)(𝑋₀−𝑥₀)^2𝛼×[𝐸_{1+𝛼,2𝛼+1}((|𝑠|^1+𝛽+|𝜆|)(𝑋₀−𝑥₀)^1+𝛼)^𝑌0_𝑦0𝑒^−𝑎0𝑡|𝑓₀(𝑡)|𝑑𝑡+(𝑋₀−𝑥₀)𝐸_{1+𝛼,2(𝛼+1)}((|𝑠|^1+𝛽+|𝜆|)(𝑋₀−𝑥₀)^1+𝛼)^𝑌0_𝑦0𝑒^−𝑎0𝑡|𝑓₁(𝑡)|𝑑𝑡], (2.19) we write solution (2.18) of the Volterra type equation (2.15) in terms of the Mittag-Leffler functions:𝑈(𝑥,𝑠)=(𝑥−𝑥₀)^𝛼−1𝐸_1+𝛼,𝛼(−(𝑠^1+𝛽−𝜆)(𝑥−𝑥₀)^1+𝛼)𝐹₀(𝑠)+(𝑥−𝑥₀)^𝛼𝐸_{1+𝛼,1+𝛼}(−(𝑠^1+𝛽−𝜆)(𝑥−𝑥₀)^1+𝛼)𝐹₁(𝑠). (2.20) In order to cancel the Laplace transformation by 𝑠 in (2.20) we will appeal to its distributional form (1.6) in Zemanian's space ℒ^′(0) (see [4]), which is dual of the countable-union space of test functions ℒ(0) defined by ℒ(0)=^∞_𝜈=1ℒ_𝑎𝜈, (2.21) where {𝑎_𝜈}^∞_𝜈=1 is a sequence of real numbers 𝑎_𝜈>0, which converges monotonically to 0+ as 𝜈→∞ and each ℒ_𝑎𝜈 is a testing-function space of smooth functions 𝜑(𝑦),𝑦∈(𝑦₀,∞) and for each nonnegative integer 𝑘 it satisfies 𝜌_𝑘,𝜈(𝜑)=sup_{𝑦∈(𝑦0,∞)}𝑒^𝑎𝜈𝑦|𝜑^(𝑘)(𝑦)|<∞,𝑘∈ℕ₀. (2.22) According to [4, Chapter III] we assign ℒ_𝑎𝜈 a topology generated by the multinorm (2.22). Consequently, ℒ_𝑎𝜈 is a countably multinormed space and the kernel of the Laplace transform 𝑒^−𝑠𝑡 is a member of ℒ_𝑎𝜈 if and only if Re𝑠≥𝑎_𝜈. Taking the space ℒ(0) we have an advantage that the space of smooth functions with compact support 𝒟 is dense in ℒ(0) and the members of the dual ℒ^′(0) are distributions. Moreover, any 𝑓∈ℒ^′(0) is a right-sided Laplace-transformable generalized function via the formula (1.6) with the right half-plane Re𝑠>0 as a region of definition. Meanwhile, any analytic function 𝐹(𝑠) on the half-plane Re𝑠≥𝑎>0, which satisfies the estimate |𝐹(𝑠)|≤𝑒^{−𝑇Re𝑠}𝑃(|𝑠|),Re𝑠≥𝑎, (2.23) where 𝑃(𝑧) is a polynomial, may be identified as the Laplace transform (1.6) of a right-sided Laplace-transformable generalized function which is concentrated on 𝑇≤𝑡<∞. Finally, the uniqueness and inversion properties are true and the inversion formula has the form 𝑓(𝑡)=lim_𝑟→∞^{𝜎+𝑖𝑟}_{𝜎−𝑖𝑟}𝐹(𝑠)𝑒^𝑠𝑡𝑑𝑠, (2.24) in the sense of convergence in 𝒟^′ for any 𝜎>0.
So in order to find eigenfunctions and general fundamental solutions of the fractional Laplace equation (2.4) we will invert the Laplace transform in (2.20) by using formula (2.24). Of course, we understand that the conventional right-sided Laplace transform (1.5) is a particular case of (1.6) being applied to a regular generalized function 𝑓∈𝐿₁((𝑇_𝑓,∞);𝑒^−𝑎0𝑡𝑑𝑡), 𝑎₀>0.
Further, we have
lim_𝑟→∞^{𝜎+𝑖𝑟}_{𝜎−𝑖𝑟}𝑠^𝑛(1+𝛽)𝐹_𝑖(𝑠)𝑒^𝑠𝑦𝑑𝑠=lim_𝑟→∞(𝑑)𝑑𝑦^{[𝑛(1+𝛽)]+1𝜎+𝑖𝑟}_{𝜎−𝑖𝑟}𝑠^{{𝑛(1+𝛽)}−1}𝐹_𝑖(𝑠)𝑒^𝑠𝑦𝑑𝑠=(𝐷_𝑦^𝑛(1+𝛽)₀𝑓_𝑖)(𝑦),𝑖=0,1,𝑛∈ℕ₀, (2.25) where {} is a fractional part of the number, the convergence is in 𝒟^′, and we assume that (𝐼_𝑦^{1−{𝑛(1+𝛽)}}₀𝑓_𝑖)(𝑦)∈𝐴𝐶^{[𝑛(1+𝛽)]+1}[𝑦₀,𝑌₀] for any 𝑛∈ℕ₀. Therefore, canceling the Laplace transformation in (2.20) and taking into account (2.16) after straightforward calculations we get the expression for a family of eigenfunctions of (2.4): 𝑢_𝜆(𝑥,𝑦)=^∞_𝑛=0(−1)^𝑛𝑛!(𝑥−𝑥₀)^{𝑛(1+𝛼)+𝛼−1}×₁Ψ₁((1,𝑛+1);(1+𝛼,𝑛(1+𝛼)+𝛼);𝜆(𝑥−𝑥₀)^1+𝛼)(𝐷_𝑦^𝑛(1+𝛽)₀𝑓₀+)(𝑦)_∞𝑛=0(−1)^𝑛𝑛!(𝑥−𝑥₀)₁^{𝑛(1+𝛼)+𝛼}Ψ₁((1,𝑛+1);(1+𝛼,(𝑛+1)(1+𝛼));𝜆(𝑥−𝑥₀)^1+𝛼)×(𝐷_𝑦^𝑛(1+𝛽)₀𝑓₁)(𝑦), (2.26) where ₁Ψ₁((𝑎₁,b₁);(𝑐₁,𝑑₁);𝑧) is the generalized Wright function [2]: ₁Ψ₁((𝑎₁,𝑏₁);(𝑐₁,𝑑₁);𝑧)=^∞_𝑚=0=Γ(𝑎₁𝑚+𝑏₁)Γ(𝑐₁𝑚+𝑑₁)𝑧^𝑚𝑚!, (2.27) and the convergence of series in (2.26) is in 𝒟^′. Letting in (2.26) 𝜆=0 we immediately come out with a classical fundamental solution of (2.4): 𝑢_𝑓(𝑥,𝑦)=^∞_𝑛=0(−1)^𝑛(𝑥−𝑥₀)^{𝑛(1+𝛼)+𝛼−1}Γ(𝑛(1+𝛼)+𝛼)(𝐷_𝑦^𝑛(1+𝛽)₀𝑓₀+)(𝑦)_∞𝑛=0(−1)^𝑛(𝑥−𝑥₀)^{𝑛(1+𝛼)+𝛼}Γ((𝑛+1)(1+𝛼))(𝐷_𝑦^𝑛(1+𝛽)₀𝑓₁)(𝑦). (2.28) Taking into account definition (1.7) of the Mittag-Leffler function, solution (2.28) may be written in the operational form 𝑢_𝑓(𝑥,𝑦)=(𝑥−𝑥₀)^𝛼−1𝐸_1+𝛼,𝛼(−(𝑥−𝑥₀)^1+𝛼𝐷_𝑦^1+𝛽₀)𝑓₀(𝑦)+(𝑥−𝑥₀)^𝛼𝐸_{1+𝛼,1+𝛼}(−(𝑥−𝑥₀)^1+𝛼𝐷_𝑦^1+𝛽₀)𝑓₁(𝑦). (2.29) Analogously, in the case of (2.10) we show that 𝑢_𝜆(𝑥,𝑦)=^∞_𝑛=0(−1)^𝑛𝑛!(𝑦−𝑦₀)^{𝑛(1+𝛽)+𝛽−1}×₁Ψ₁((1,𝑛+1);(1+𝛽,𝑛(1+𝛽)+𝛽);𝜆(𝑦−𝑦₀)^1+𝛽)(𝐷_𝑥^𝑛(1+𝛼)₀ℎ₀+)(𝑥)_∞𝑛=0(−1)^𝑛𝑛!(𝑦−𝑦₀)₁^{𝑛(1+𝛽)+𝛽}Ψ₁((1,𝑛+1);(1+𝛽,(𝑛+1)(1+𝛽));𝜆(𝑦−𝑦₀)^1+𝛽)×(𝐷_ℎ^𝑛(1+𝛼)₀ℎ₁𝑢)(𝑥),_ℎ(𝑥,𝑦)=(𝑦−𝑦₀)^𝛽−1𝐸_1+𝛽,𝛽(−(𝑦−𝑦₀)^1+𝛽𝐷_𝑥^1+𝛼₀)ℎ₀(𝑥)+(𝑦−𝑦₀)^𝛽𝐸_{1+𝛽,1+𝛽}(−(𝑦−𝑦₀)^1+𝛽𝐷_𝑥^1+𝛼₀)ℎ₁(𝑥) (2.30)(2.31) are also correspondingly eigenfunctions and a fundamental solution of (2.4).
On the other hand we may write solutions (2.29) and (2.31) in the form of the generalized Neumann series. Namely, we find
𝑢_𝑓1(𝑥,𝑦)=Γ(𝛼)^∞_𝑛=0(−1)^𝑛𝐼_𝑥^𝑛(1+𝛼)₀(𝑥−𝑥₀)^𝛼−1(𝐷_𝑦^𝑛(1+𝛽)₀𝑓₀+1)(𝑦)Γ(1+𝛼)^∞_𝑛=0(−1)^𝑛𝐼_𝑥^𝑛(1+𝛼)₀(𝑥−𝑥₀)^𝛼(𝐷_𝑦^𝑛(1+𝛽)₀𝑓₁)(𝑦), (2.32) taking in mind the analyticity of series in (2.32) by 𝑥 in the interval 0<𝜀≤|𝑥−𝑥₀|<𝑅 and by 𝑦∈(𝑦₀,𝑌₀). In the same manner, we represent (2.31) by the expression 𝑢_ℎ1(𝑥,𝑦)=Γ(𝛽)^∞_𝑛=0(−1)^𝑛𝐼_𝑦^𝑛(1+𝛽)₀(𝑦−𝑦₀)^𝛽−1(𝐷_𝑥^𝑛(1+𝛼)₀ℎ₀+1)(𝑥)Γ(1+𝛽)^∞_𝑛=0(−1)^𝑛𝐼_𝑦^𝑛(1+𝛽)₀(𝑦−𝑦₀)^𝛽(𝐷_𝑥^𝑛(1+𝛼)₀ℎ₁)(𝑥) (2.33) with arbitrary ℎ₀,ℎ₁ assuming analyticity of the corresponding series by 𝑦 in the interval 0<𝜀≤|𝑦−𝑦₀|<𝑅 and by 𝑥∈(𝑥₀,𝑋₀). Now taking into account zero values 𝐷_𝑥^1+𝛼₀₊(𝑥−𝑥₀)^𝛼−1=0, 𝐷_𝑥^1+𝛼₀₊(𝑥−𝑥₀)^𝛼=0, 𝐷_𝑦^1+𝛽₀₊(𝑦−𝑦₀)^𝛽−1=0, 𝐷_𝑦^1+𝛽₀₊(𝑦−𝑦₀)^𝛼=0 it is not difficult to verify that (2.32) and (2.33) are classical fundamental solutions of the fractional Laplacian Δ^𝛼,𝛽𝑢=0 subject to conditions (2.11), (2.12), (2.13), and (2.14) respectively. Thus we have proved

Theorem 2.6. In case (ii) functions (2.26) and (2.30) represent eigenfunctions of the fractional Laplacian (2.4) and expressions (2.28) and (2.31) are unique classical fundamental solutions subject to conditions (2.11), (2.12), (2.13), and (2.14) respectively. These solutions can be written in the corresponding form of generalized Neumann series (2.32) and (2.33) under additional conditions of analyticity.

Finally, in case (iii) an analog of (2.9) and (2.10) is

(𝐼_𝑥^1+𝛼₀₊𝑢)(𝑥,𝑦)+(𝐼_𝑦^1+𝛽₀₊𝑢)(𝑥,𝑦)−𝜆(𝐼_𝑥^1+𝛼₀₊𝐼_𝑦^1+𝛽₀₊=𝑢)(𝑥,𝑦)(𝑥−𝑥₀)^𝛼−1Γ(𝛼)(𝐼_𝑦^1+𝛽₀₊𝑓₀)(𝑦)+(𝑥−𝑥₀)^𝛼Γ(1+𝛼)(𝐼_𝑦^1+𝛽₀₊𝑓₁+)(𝑦)(𝑦−𝑦₀)^𝛽−1Γ(𝛽)(𝐼_𝑥^1+𝛼₀₊ℎ₀)(𝑥)+(𝑦−𝑦₀)^𝛽Γ(1+𝛽)(𝐼_𝑥^1+𝛼₀₊ℎ₁)(𝑥). (2.34) Consequently, in the right-hand side of Volterra's equation (2.15) we get an additional term

𝑈(𝑥,𝑠)+(𝑠^1+𝛽−𝜆)Γ(1+𝛼)^𝑥_𝑥0(𝑥−𝑡)^𝛼𝑈(𝑡,𝑠)𝑑𝑡=𝐹(𝑥,𝑠)+𝑒^−𝑦0𝑠(𝐼_𝑥^1+𝛼₀₊(𝑠ℎ₀+ℎ₁))(𝑥), (2.35) which will give a source for generalized eigenfunctions and fundamental solutions of the fractional Laplacian (2.4). In fact, owing to the estimates

|(𝑠^1+𝛽−𝜆)^𝑥_𝑥0(𝑥−𝑡)^𝛼𝐸_{1+𝛼,1+𝛼}(−(𝑠^1+𝛽−𝜆)(𝑥−𝑡)^1+𝛼)𝐹(𝑡,𝑠)𝑑𝑡|≤|𝐹₀(𝑠)|^∞_𝑛=0(|𝑠|^1+𝛽+|𝜆|)^𝑛+1Γ((𝑛+1)(1+𝛼))Γ(𝛼)^𝑥_𝑥0(𝑥−𝑡)^{𝑛(1+𝛼)+𝛼}(𝑡−𝑥₀)^𝛼−1𝑑𝑡+|𝐹₁(𝑠)|^∞_𝑛=0(|𝑠|^1+𝛽+|𝜆|)^𝑛+1Γ((𝑛+1)(1+𝛼))Γ(1+𝛼)^𝑥_𝑥0(𝑥−𝑡)^{𝑛(1+𝛼)+𝛼}(𝑡−𝑥₀)^𝛼+𝑑𝑡_∞𝑛=0|𝑠|(|𝑠|^1+𝛽+|𝜆|)^𝑛+1Γ((𝑛+1)(1+𝛼))Γ(1+𝛼)^𝑥_𝑥0(𝑥−𝑡)^{𝑛(1+𝛼)+𝛼𝑡}_𝑥0(𝑡−𝑣)^𝛼|ℎ₀+(𝑣)|𝑑𝑣𝑑𝑡_∞𝑛=0(|𝑠|^1+𝛽+|𝜆|)^𝑛+1Γ((𝑛+1)(1+𝛼))Γ(1+𝛼)^𝑥_𝑥0(𝑥−𝑡)^{𝑛(1+𝛼)+𝛼𝑡}_𝑥0(𝑡−𝑣)^𝛼|ℎ₁(𝑣)|𝑑𝑣𝑑𝑡≤(|𝑠|^1+𝛽+|𝜆|)(𝑋₀−𝑥₀)^2𝛼[𝐸_{1+𝛼,2𝛼+1}((|𝑠|^1+𝛽+|𝜆|)(𝑋₀−𝑥₀)^1+𝛼)^𝑌0_𝑦0𝑒^−𝑎0𝑡|𝑓₀(𝑡)|𝑑𝑡+(𝑋₀−𝑥₀)𝐸_{1+𝛼,2(𝛼+1)}((|𝑠|^1+𝛽+|𝜆|)(𝑋₀−𝑥₀)^1+𝛼)×[^𝑌0_𝑦0𝑒^−𝑎0𝑡|𝑓₁(𝑡)|𝑑𝑡+|𝑠|^𝑋0_𝑥0|ℎ₀(𝑡)|𝑑𝑡+^𝑋0_𝑥0|ℎ₁(𝑡)|𝑑𝑡]]<∞, (2.36) we write solution of the Volterra type equation (2.35) in terms of the Mittag-Leffler functions and generalized Neumann series:

𝑈(𝑥,𝑠)=(𝑥−𝑥₀)^𝛼−1𝐸_1+𝛼,𝛼(−(𝑠^1+𝛽−𝜆)(𝑥−𝑥₀)^1+𝛼)𝐹₀(𝑠)+(𝑥−𝑥₀)^𝛼𝐸_{1+𝛼,1+𝛼}(−(𝑠^1+𝛽−𝜆)(𝑥−𝑥₀)^1+𝛼)𝐹₁(𝑠)+𝑒^{−𝑦0𝑠∞}_𝑛=0(−1)^𝑛(𝑠^1+𝛽−𝜆)^𝑛(𝐼_𝑥^{(𝑛+1)(1+𝛼)}₀(𝑠ℎ₀+ℎ₁))(𝑥). (2.37) Cancelling the Laplace transformation we take in mind the relations

lim_𝑟→∞^{𝜎+𝑖𝑟}_{𝜎−𝑖𝑟}𝑠^𝑛(1+𝛽)𝑒^{𝑠(𝑦−𝑦0)}𝑑𝑠=lim_𝑟→∞(𝑑)𝑑𝑦^{[𝑛(1+𝛽)]+1𝜎+𝑖𝑟}_{𝜎−𝑖𝑟}⟨𝑡₊^{−{𝑛(1+𝛽)}}Γ(1−{𝑛(1+𝛽)})∗𝛿(𝑡−𝑦₀),𝑒^−𝑠𝑡⟩𝑒^𝑠𝑦𝑑𝑑𝑠=()𝑑𝑦^{[𝑛(1+𝛽)]+1}(𝑡₊^{−{𝑛(1+𝛽)}}Γ(1−{𝑛(1+𝛽)})∗𝛿(𝑡−𝑦₀))(𝑦)=(𝐷_𝑦^𝑛(1+𝛽)₀𝛿)(𝑦), (2.38) where 𝛿 is Dirac's delta-function and ∗ denotes the convolution product. Therefore, after straightforward calculations we get the expression for a family of eigenfunctions of (2.4):

𝑢_𝜆(𝑥,𝑦)=^∞_𝑛=0(−1)^𝑛𝑛!(𝑥−𝑥₀)^{𝑛(1+𝛼)+𝛼−1}×₁Ψ₁((1,𝑛+1);(1+𝛼,𝑛(1+𝛼)+𝛼);𝜆(𝑥−𝑥₀)^1+𝛼)(𝐷_𝑦^𝑛(1+𝛽)₀𝑓₀+)(𝑦)_∞𝑛=0(−1)^𝑛𝑛!(𝑥−𝑥₀)₁^{𝑛(1+𝛼)+𝛼}Ψ₁((1,𝑛+1);(1+𝛼,(𝑛+1)(1+𝛼));𝜆(𝑥−𝑥₀)^1+𝛼)×(𝐷_𝑦^𝑛(1+𝛽)₀𝑓₁)(𝑦)+^𝑥_𝑥0ℎ₁(𝑡)^∞_𝑛=0(−1)^𝑛𝑛!×(𝑥−𝑡)₁^{𝑛(1+𝛼)+𝛼}Ψ₁((1,𝑛+1);(1+𝛼,(𝑛+1)(1+𝛼));𝜆(𝑥−𝑡)^1+𝛼)×(𝐷_𝑦^𝑛(1+𝛽)₀𝛿)(𝑦)𝑑𝑡+^𝑥_𝑥0ℎ₀(𝑡)^∞_𝑛=0(−1)^𝑛𝑛!(𝑥−𝑡)^{𝑛(1+𝛼)+𝛼}×₁Ψ₁((1,𝑛+1);(1+𝛼,(𝑛+1)(1+𝛼));𝜆(𝑥−𝑡)^1+𝛼)×(𝐷_𝑦^𝑛(1+𝛽)₀𝛿^′)(𝑦)𝑑𝑡, (2.39) where the convergence of series in (2.39) is in 𝒟^′. Letting in (2.39) 𝜆=0 we derive a generalized fundamental solution of (2.4):

𝑢(𝑥,𝑦)=^∞_𝑛=0(−1)^𝑛(𝑥−𝑥₀)^{𝑛(1+𝛼)+𝛼−1}Γ(𝑛(1+𝛼)+𝛼)(𝐷_𝑦^𝑛(1+𝛽)₀𝑓₀+)(𝑦)_∞𝑛=0(−1)^𝑛(𝑥−𝑥₀)^{𝑛(1+𝛼)+𝛼}Γ((𝑛+1)(1+𝛼))(𝐷_𝑦^𝑛(1+𝛽)₀𝑓₁)(𝑦)+^𝑥_𝑥0ℎ₁(𝑡)^∞_𝑛=0(−1)^𝑛(𝑥−𝑡)^{(𝑛+1)(1+𝛼)−1}Γ((𝑛+1)(1+𝛼))(𝐷_𝑦^𝑛(1+𝛽)₀+𝛿)(𝑦)𝑑𝑡_𝑥𝑥0ℎ₀(𝑡)^∞_𝑛=0(−1)^𝑛(𝑥−𝑡)^{(𝑛+1)(1+𝛼)−1}Γ((𝑛+1)(1+𝛼))(𝐷_𝑦^𝑛(1+𝛽)₀𝛿^′)(𝑦)𝑑𝑡, (2.40) which may be written in the operational form

𝑢(𝑥,𝑦)=(𝑥−𝑥₀)^𝛼−1𝐸_1+𝛼,𝛼(−(𝑥−𝑥₀)^1+𝛼𝐷_𝑦^1+𝛽₀)𝑓₀(𝑦)+(𝑥−𝑥₀)^𝛼𝐸_{1+𝛼,1+𝛼}(−(𝑥−𝑥₀)^1+𝛼𝐷_𝑦^1+𝛽₀)𝑓₁+(𝑦)_𝑥𝑥0(𝑥−𝑡)^𝛼𝐸_{1+𝛼,1+𝛼}(−(𝑥−𝑡)^1+𝛼𝐷_𝑦^1+𝛽₀)𝛿(𝑦)ℎ₁(𝑡)𝑑𝑡+^𝑥_𝑥0(𝑥−𝑡)^𝛼𝐸_{1+𝛼,1+𝛼}(−(𝑥−𝑡)^1+𝛼𝐷_𝑦^1+𝛽₀)𝛿^′(𝑦)ℎ₀(𝑡)𝑑𝑡. (2.41) Analogously, functions

𝑢_𝜆(𝑥,𝑦)=^∞_𝑛=0(−1)^𝑛𝑛!(𝑦−𝑦₀)^{𝑛(1+𝛽)+𝛽−1}×₁Ψ₁((1,𝑛+1);(1+𝛽,𝑛(1+𝛽)+𝛽);𝜆(𝑦−𝑦₀)^1+𝛽)(𝐷_𝑥^𝑛(1+𝛼)₀ℎ₀+)(𝑥)_∞𝑛=0(−1)^𝑛𝑛!(𝑦−𝑦₀)₁^{𝑛(1+𝛽)+𝛽}Ψ₁((1,𝑛+1);(1+𝛽,(𝑛+1)(1+𝛽));𝜆(𝑦−𝑦₀)^1+𝛽)×(𝐷_ℎ^𝑛(1+𝛼)₀ℎ₁)(𝑥)+^𝑦_𝑦0ℎ₁(𝑡)^∞_𝑛=0(−1)^𝑛𝑛!×(𝑦−𝑡)₁^{𝑛(1+𝛽)+𝛽}Ψ₁((1,𝑛+1);(1+𝛽,(𝑛+1)(1+𝛽));𝜆(𝑦−𝑡)^1+𝛽)×(𝐷_𝑥^𝑛(1+𝛼)₀𝛿)(𝑥)𝑑𝑡+^𝑦_𝑦0ℎ₀(𝑡)^∞_𝑛=0(−1)^𝑛𝑛!(𝑦−𝑡)^{𝑛(1+𝛽)+𝛽}×₁Ψ₁((1,𝑛+1);(1+𝛽,(𝑛+1)(1+𝛽));𝜆(𝑦−𝑡)^1+𝛽)×(𝐷_𝑥^𝑛(1+𝛼)₀𝛿^′)(𝑥)𝑑𝑡,𝑢(𝑥,𝑦)=(𝑦−𝑦₀)^𝛽−1𝐸_1+𝛽,𝛽(−(𝑦−𝑦₀)^1+𝛽𝐷_𝑥^1+𝛼₀)ℎ₀(𝑥)+(𝑦−𝑦₀)^𝛽𝐸_{1+𝛽,1+𝛽}(−(𝑦−𝑦₀)^1+𝛽𝐷_𝑥^1+𝛼₀)ℎ₁+(𝑥)_𝑦𝑦0(𝑦−𝑡)^𝛽𝐸_{1+𝛽,1+𝛽}(−(𝑦−𝑡)^1+𝛽𝐷_𝑥^1+𝛼₀)𝛿(𝑥)𝑓₁(𝑡)𝑑𝑡+^𝑦_𝑦0(𝑦−𝑡)^𝛽𝐸_{1+𝛽,1+𝛽}(−(𝑦−𝑡)^1+𝛽𝐷_𝑥^1+𝛼₀)𝛿^′(𝑥)𝑓₀(𝑡)𝑑𝑡 (2.42)(2.43) are also correspondingly eigenfunctions and generalized fundamental solutions of (2.4).

Theorem 2.7. In case (iii) functions (2.39) and (2.42) represent eigenfunctions and expressions (2.40), (2.41), and (2.43) are generalized fundamental solutions of the fractional Laplacian (2.4).

Example 2.8. As a particular case, it is not difficult to obtain from (2.28), (2.31) the classical fundamental solution 𝑢(𝑥,𝑦)=(1/2)log((𝑥−𝑥₀)²+(𝑦−𝑦₀)²) of the Laplace equation Δ𝑢=0(𝛼=𝛽=1). Indeed, putting 𝑓₁(𝑦)=ℎ₁(𝑥)=0, 𝑓₀(𝑦)=log(𝑦−𝑦₀), ℎ₀(𝑥)=log(𝑥−𝑥₀) we assume, correspondingly, (𝑥−𝑥₀)/(𝑦−𝑦₀)<1, (𝑥−𝑥₀)/(𝑦−𝑦₀)>1 in (2.28), (2.31), and for instance, solution 𝑢_𝑓(𝑥,𝑦) becomes 𝑢_𝑓(𝑥,𝑦)=^∞_𝑛=0(−1)^𝑛(𝑥−𝑥₀)^2𝑛(2𝑛)!(log(𝑦−𝑦₀))^(2𝑛)=log(𝑦−𝑦₀)+^∞_𝑛=1(−1)^𝑛+1(2𝑛𝑥−𝑥₀𝑦−𝑦₀)^2𝑛=log(𝑦−𝑦₀)+^∞_𝑛=0(−1)^𝑛(2(𝑛+1)𝑥−𝑥₀𝑦−𝑦₀)^2(𝑛+1)=log(𝑦−𝑦₀1)+2log(1+(𝑥−𝑥₀𝑦−𝑦₀)²)=12log((𝑥−𝑥₀)²+(𝑦−𝑦₀)²). (2.44) Analogously we treat solution (2.31).

3. Separation of Variables: Analytic Solutions

The method of separation of variables allows us to simplify eigenfunctions and fundamental solutions of the fractional Laplacian. Indeed, putting 𝑢(𝑥,𝑦)=𝑢₁(𝑥)𝑢₂(𝑦), substituting in (2.34), and taking into account initial conditions (2.11), (2.12), (2.13), and (2.14) it becomes

𝑢₂(𝑦)(𝐼_𝑥^1+𝛼₀₊𝑢₁)(𝑥)+𝑢₁(𝑥)(𝐼_𝑦^1+𝛽₀₊𝑢₂)(𝑦)−𝜆(𝐼_𝑥^1+𝛼₀₊𝑢₁)(𝑥)(𝐼_𝑦^1+𝛽₀₊𝑢₂)(𝑦)=𝑎₁(𝑥−𝑥₀)^𝛼−1Γ(𝛼)(𝐼_𝑦^1+𝛽₀₊𝑢₂)(𝑦)+𝑎₂(𝑥−𝑥₀)^𝛼Γ(1+𝛼)(𝐼_𝑦^1+𝛽₀₊𝑢₂)(𝑦)+𝑏₁(𝑦−𝑦₀)^𝛽−1Γ(𝛽)(𝐼_𝑥^1+𝛼₀₊𝑢₁)(𝑥)+𝑏₂(𝑦−𝑦₀)^𝛽Γ(1+𝛽)(𝐼_𝑥^1+𝛼₀₊𝑢₁)(𝑥), (3.1) where 𝑎_𝑖,𝑏_𝑖∈ℂ, 𝑖=1,2 are arbitrary constants. If (𝐼_𝑥^1+𝛼₀₊𝑢₁)(𝑥)(𝐼_𝑦^1+𝛽₀₊𝑢₂)(𝑦)≠0, (𝑥,𝑦)∈(𝑥₀,𝑋₀)×(𝑦₀,𝑌₀) we divide (3.1) by this product and separate variables, getting two Abel's type second kind integral equations to define 𝑢_𝑖, 𝑖=1,2, namely,

𝑢₁(𝑥)+𝜇(𝐼_𝑥^1+𝛼₀₊𝑢₁)(𝑥)=𝑎₁(𝑥−𝑥₀)^𝛼−1Γ(𝛼)+𝑎₂(𝑥−𝑥₀)^𝛼,𝑢Γ(1+𝛼)₂(𝑦)−(𝜆+𝜇)(𝐼_𝑦^1+𝛽₀₊𝑢₂)(𝑦)=𝑏₁(𝑦−𝑦₀)^𝛽−1Γ(𝛽)+𝑏₂(𝑦−𝑦₀)^𝛽,Γ(1+𝛽) (3.2) where 𝜆,𝜇∈ℂ are constants. We note that the equality (𝐼_𝑥^1+𝛼₀₊𝑢₁)(𝑥)(𝐼_𝑦^1+𝛽₀₊𝑢₂)(𝑦)=0 for at least one point (𝜉,𝜂) agrees with (3.1) and (3.2). So we solve the latter equations similarly to (2.18), arriving at the following family of eigenfunctions 𝑢_𝜆,𝜇(𝑥,𝑦)=𝑢₁(𝑥)𝑢₂(𝑦), where

𝑢₁(𝑥)=𝑎₁(𝑥−𝑥₀)^𝛼−1𝐸_1+𝛼,𝛼(𝜇(𝑥−𝑥₀)^1+𝛼)+𝑎₂(𝑥−𝑥₀)^𝛼𝐸_{1+𝛼,1+𝛼}(𝜇(𝑥−𝑥₀)^1+𝛼𝑢),₂(𝑦)=𝑏₁(𝑦−𝑦₀)^𝛽−1𝐸_1+𝛽,𝛽((𝜆−𝜇)(𝑦−𝑦₀)^1+𝛽)+𝑏₂(𝑦−𝑦₀)^𝛽𝐸_{1+𝛽,1+𝛽}((𝜆−𝜇)(𝑦−𝑦₀)^1+𝛽). (3.3) On the other hand, we may write these solutions in terms of the generalized Neumann series. Precisely, denoting by

𝑣_𝛼(𝑥)=𝑎₁(𝑥−𝑥₀)^𝛼−1Γ(𝛼)+𝑎₂(𝑥−𝑥₀)^𝛼,𝑣Γ(1+𝛼)_𝛽(𝑦)=𝑏₁(𝑦−𝑦₀)^𝛽−1Γ(𝛽)+𝑏₂(𝑦−𝑦₀)^𝛽,Γ(1+𝛽) (3.4) and recalling index properties for the fractional integral (1.2) we get representations of (3.3) in the respective resolvent form for fractional integral operators 𝐼_𝑥^1+𝛼₀∶𝐿₁(𝑥₀,𝑋₀)→𝐿₁(𝑥₀,𝑋₀), 𝐼_𝑦^1+𝛽₀∶𝐿₁(𝑦₀,𝑌₀)→𝐿₁(𝑦₀,𝑌₀):

𝑢₁(𝑥)=^∞_𝑛=0𝐼_𝑥^𝑛(1+𝛼)₀𝑣_𝛼𝜇^𝑛=(𝐸−𝜇𝐼_𝑥^1+𝛼₀)⁻¹𝑣_𝛼=ℛ(𝜇;𝐼_𝑥^1+𝛼₀)𝑣_𝛼,𝑢₂(𝑦)=^∞_𝑛=0𝐼_𝑦^𝑛(1+𝛽)₀𝑣_𝛽(𝜆−𝜇)^𝑛=(𝐸−(𝜆−𝜇)𝐼_𝑦^1+𝛽₀)⁻¹𝑣_𝛽=ℛ(𝜆−𝜇;𝐼_𝑦^1+𝛽₀)𝑣_𝛽, (3.5) where 𝐸 as usual is the identity operator. It is easily seen that series (3.5) are analytic with respect to 𝑥∈(𝑥₀,𝑋₀) and 𝑦∈(𝑦₀,𝑌₀). Further, since (see (1.2))

‖𝐼_𝑥^1+𝛼₀𝑓‖_{𝐿1(𝑥0,𝑋0)}≤(𝑋₀−𝑥₀)^𝛼+1Γ(1+𝛼)‖𝑓‖_{𝐿1(𝑥0,𝑋0)}, (3.6) we have

‖ℛ(𝑧;𝐼_𝑥^1+𝛼₀)‖≤^∞_𝑛=0‖𝐼_𝑥^1+𝛼₀‖^𝑛|𝑧|^𝑛≤^∞_𝑛=0((𝑋₀−𝑥₀)^𝛼+1)Γ(1+𝛼)^𝑛|𝑧|^𝑛. (3.7) Therefore the resolvent functions ℛ(𝑧;𝐼_𝑥^1+𝛼₀), ℛ(𝑧;𝐼_𝑦^1+𝛽₀) are analytic in the open discs:

|𝑧|<Γ(1+𝛼)(𝑋₀−𝑥₀)^𝛼+1,|𝑧|<Γ(1+𝛽)(𝑌₀−𝑦₀)^𝛽+1, (3.8) respectively. Thus we write a family of eigenfunctions for (2.4) in the resolvent form

𝑢_𝜆(𝑥,𝑦)=ℛ(𝜇;𝐼_𝑥^1+𝛼₀)𝑣_𝛼ℛ(𝜆−𝜇;𝐼_𝑦^1+𝛽₀)𝑣_𝛽. (3.9) Indeed, substituting (3.9) into (2.4) taking into account the values 𝐷_𝑥^1+𝛼₀₊𝑣_𝛼=𝐷_𝑦^1+𝛽₀₊𝑣_𝛽=0 after a simple change of the summation index into the series we easily satisfy (2.4). But we will extend our family of eigenfunctions considering

𝑢_𝜆(𝑥,𝑦)=ℛ(𝜇;𝐼_𝑥^1+𝛼₀)𝑓ℛ(𝜆−𝜇;𝐼_𝑦^1+𝛽₀)𝑔,𝜇∈ℂ, (3.10) with arbitrary 𝑓(𝑥),𝑔(𝑦) such that 𝐼_𝑥^{𝑛(1+𝛼)+1−𝛼}₀𝑓∈𝐴𝐶²([𝑥₀,𝑋₀]), 𝐼_𝑦^{𝑛(1+𝛽)+1−𝛽}₀𝑔∈𝐴𝐶²([𝑦₀,𝑌₀]), and the corresponding resolvent function (3.5) are analytic by 𝑥 and 𝑦. So substituting (3.10) into (2.4) and ignoring trivial cases 𝐷_𝑥^1+𝛼₀₊𝑓=𝐷_𝑦^1+𝛽₀₊𝑔=0 which drive immediately to (3.10) with 𝑓=𝑣_𝛼, 𝑔=𝑣_𝛽 (see (2.3)), after separation of variables we obtain fractional differential equations to define 𝑓,𝑔:

𝐷_𝑥^1+𝛼₀₊𝑓=𝑐ℛ(𝜇;𝐼_𝑥^1+𝛼₀𝐷)𝑓,_𝑦1+𝛽0+𝑔=−𝑐ℛ(𝜆−𝜇;𝐼_𝑦^1+𝛽₀)𝑔, (3.11) where 𝑐 is an arbitrary constant. Hence acting by inverse operators 𝐸−𝜇𝐼_𝑥^1+𝛼₀ and 𝐸−(𝜆−𝜇)𝐼_𝑦^1+𝛽₀ on (3.11) with the use of (2.3) we get, correspondingly,

𝐷_𝑥^1+𝛼₀₊𝑓=(𝜇+𝑐)𝑓(𝑥)+𝜇𝑣_𝛼,𝐷_𝑦^1+𝛽₀₊𝑔=(𝜆−𝜇−𝑐)𝑔(𝑦)+(𝜆−𝜇)𝑣_𝛽. (3.12) The latter equations are solved, for instance, in [1, 2] and we obtain the following solutions (𝑐_𝑖,𝑑_𝑖,𝑖=1,2 are constants):

𝑓(𝑥)=𝑐₁(𝑥−𝑥₀)^𝛼−1𝐸_1+𝛼,𝛼((𝜇+𝑐)(𝑥−𝑥₀)^1+𝛼)+𝑐₂(𝑥−𝑥